ADA050879 


m /dso/.4-/i^ 


REPORT;  J 

mJMERICAL  ^DIES  OF  SEVERAL  ^DAPTIVE 

ITERATIVE  ALGORITHMs’t 

■ ■ - * * 


by 


iS-  David  R.^Kincaid 

j.T)Aug~»«77] 


ITR-CNA-126 


] 


J-S  \<P fi H C.Q>^l -'! i 


Abstract 


Six  adaptive  iterative  algorithms  are  studied  for  six  elliptic  partial 
differential  equations  on  six  regions  compatible  with  subroutine  REGION. 
An  effort  was  made  to  make  the  resulting  preliminary  ITPACK  code  conform 
to  the  "ELLPACK  Contributor's  Guide--Initial  Version^"  CSD  TR  208, 

Purdue  University,  November  1,  1976. 
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1 .  Int  roduc  t ion 

The  initial  ITPACK  code  was  conceived  to  be  a package  of 
Fortran  subroutines  to  solve  the  large  sparse  positive  definite 
linear  systems  which  arise  from  the  five-point  finite  difference 
Jisoreti aacion  of  a general  self-adjoint  elliptic  partial  differen- 
cial equation 

(1.1)  (au)  +(cu)  +fu=g 

X X y y ** 

with  Dirichlet  boundary  conditions  on  a region  compatible  with  the 
REGION  subprogram.  (For  details,  see  Appendix  3.) 

The  current  ITPACK  code  contains  the  following  six  iterative 
algorithms 

I.  Jacobi  Semi-iteration  (J-SI) 

II.  Compressed  Jacobi  Conjugate  Gradient  (CJ-CG) 

III.  Reduced  System  Semi-iteration  (RS-SI) 

IV.  Reduced  System  Conjugate  Gradient  (RS-CG) 

V.  Symmetric  Successive  Overrelaxation  Semi- iterat ion 
(SSOR-SI) 

VI.  Symmetric  Successive  Overrelaxation  Conjugate  Gradient 
(SSOR-SI) 

which  are  developed  in  the  monograph  [1]  by  Hageman  and  Young. 

These  methods  will  not  be  motivated  in  this  report;  however,  detailed 
algorithms  are  given  in  Appendix  1. 

Future  plans  for  the  ITPACK  project  are  many  and  varied  with 
the  major  limiting  factor  being  time  for  implementation  of  the  code. 
Various  other  iterative  algorithms  are  being  considered  at  this  time. 
These  include  the  Block  Jacobi  Semi-iterative  Method  and  the  Block 
Conjugate  Gradient  Method.  Also  coding  schemes  for  mixed  and  Neumann 
boundary  conditions  are  being  developed.  Yet  another  phase  of  this 
project  is  the  use  of  various  finite  difference  stencils . 

The  purpose  of  the  ITPACK  project  is  to  develop,  study,  and 
analyze  iterative  algoirthms  for  solving  elliptic  partial  differ- 
ential equations.  The  principal  activities  are  centered  around 
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iraproving  those  iterative  algorithms  w hich  involve  efficient 
stopping  tests  and  effective  parameter  determination  when  computing 
the  numerical  solution  of  the  large  sparse  matrix  problems  from 
elliptic  equations  whenever  finite  difference  or  finite  element 
procedures  are  employed. 

It  is  anticipated  that  the  code  from  the  ITPACK  project  will 


...  ve 


following  benefits  and  utilization: 


(a) 

(b) 

(c) 

(d) 


the  development  of  ELLPACK  modules  which  use  adaptive 
iterative  procedures  to  solve  the  linear  systems 
add  to  existing  knowledge  of  the  effectiveness  of  various 
iterative  algorithms 

allow  comparisons  between  these  iterative  schemes  and 
between  iterative  and  direct  methods 

the  development  of  quality  software  as  a research  and 
teaching  tool 

In  Section  2,  background  material  relating  (1.1)  and  basic 
iterative  methods  is  given  with  the  detailed  adaptive  iterative 
algorithms  stated  in  Appendix  1.  The  overall  structure  of  ITPACK 
is  outlined  in  Section  3 with  a sample  of  the  code  required  to 
use  the  current  ITPACK.  The  six  test  problems  and  "St  regions 

are  set-forth  in  Section  4 with  complete  details  on  subroutine 
REGION  in  Appendix  3.  Numerical  results  and  figures  are  given  in 
Section  5. 


2 . Iterative  Methods 

The  six  iterative  algorithms  covered  by  this  study  are 
developed  in  the  monograph  [ 1]  by  Hageman  and  Young.  Consequently, 
we  will  not  repeat  these  derivations  here.  We  will,  however,  pre- 
sent some  material  related  to  the  development  of  these  algorithms 
which  will  aid  in  the  understanding  of  the  detailed  statements  of 
those  procedures  given  in  Appendix  1. 

We  consider  the  general  self-adjoint  elliptic  partial  differen- 
tial aquation  with  Dirich let  boundary  conditions. 


J 
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(2.1) 


X X 


y y 


fu 

= g 

, (x,y)  e R 

u 

“ q 

, (x,y)  c 3r 

of 

both 

X and  y,  and  the  region 

is  denoced  R with  boundary  9R. 

Using  the  central  difference  discretization  at  the  grid  point 
associated  with  (i,j),  we  have 


« { (au^  ) 


“ (au_) 


(i.j)  ""  (i+Ss,j)  * (i-^.j) 


}h 


-1 


-1 


Hence,  we  use 


<“«>x  "l+l.i  "i-l.j 


-2 


Similarly , 
(cu 


we 


y)y 


use 

I 

(i 


j+ia“i,l+l  ‘^i,  j-Js“i.  j-1 


Here  we  let  u^^  denote  the  discrete  variable  as  opposed  to  the 
continuous  variable  u. 

'I'll  II  a at  (i>J),  thu  nui  l -<ul.)i>  Lut  uiilptlc  e<iuuClou  (2.1)  lb 
a P |.t  L'o  x f Ilia  c lul  by  tliu  line  nr  eqiintlon 


(2.2) 


‘“ij^'i+l.j  ‘■''ij^l , j-t-1 
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From  (2.3),  we  have  the  following  symmetry  condition 


W.  . = E.  , . 

i-l.J 


S . . = N . . , 

13  1.3-1 


so  that  only  f our coef f ic ieit  values  need  to  be  stored  per  grid  point, 
Hence,  we  have 


(2.4) 


- E.  u.,,  . - N..U.  = R.. 

13  1+1,3  13  1.3+1  13 


@ (i.3 ) 


'i-l.j 


i.3-1 


Since  only  regular  grid  points  are  considered,  we  have  for 
the  basic  linear  equation 


'^ij  ■ ^^ij'^i+l.j  ^i3^i.J  + l ^i-l.j'^i-l.j 


(2.5) 


+ N.  . ,u,  . + R,.)/C,, 

i.J-1  l.J-1  Ij  IJ 
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to 


and 


Using  matrix  notation,  equations  (2.4)  and  (2.5)  correspond 


Au  = b 


(2.6)  u = Bu  + c 


.-1 


respectively,  where  D "A  = I-B  and  D=diag(C^^).  Notice  that  if 
the  k-th  equation  in  Au  = b corresponds  to  the  grid-point  (i,j) 
then  b,  is  equal  to  R..  plus  the  sum  of  some  terms  in  (2.2), 

X 1 j ’ 

with  u replaced  by  q^  for  boundary-points  adjacent  to  (i,j). 
Clearly,  A is  symmetric  while  B is  not.  It  can  be  shown  that 
A is  positive  definite. 

When  the  red-black  ordering  is  used,  the  basic  iterative 
system  (2.6)  assumes  the  form 


'R 


'B  I 


(2.8) 


(n+1) 
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or  ia  matrix  form 


(2.9) 


+ c 


For  the  reduced  system,  the  basic  iterative  equation  would  be 


ufn+l) 


(E..uf^t^^^  + N.  + E + N 

ij  i,j+l  i“l,j  i-l,j  * i , j - l^i  , j ^ ]. 


(2.10) 


where 


(n+h) 


kA'"k+i,Jl  ^\ji'"k,i+l  ■*■  ^k-l,Jl“k-l,Ji  ”k,£-l^k?i-l 


This  corresponds  to  the  following  stencil  at  each  regular  Iterior 
grid  point 


t'‘i.j+2 

u 


i+1, j+1 


u . _ . 

1-2.  J 


i-1, j-1 


''i+2,j 


L.  ,, 

1.  J-2 


""i+l.j-l 


The  basic  iterative  equation  for  the  resulted  system  is 


(2.11)  - FjF^u<”>  + t Cj 


However,  it  is  easier  to  consider  this  as  two  separate  iterations. 
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(n+ii)  ^ „ (n)  ^ ^ 


= F + 


First,  ii  sweep  of  the  red  grid  points  would  be  done  which  involves 
a "weighted-average"  of  the  adjacent  black  grid  points  with  the 
results  being  stored  in  the  red  storage  locations.  This  is  followed 
by  a similar  sweep  of  the  black  grid  points  using  (2.10).  The  net 
result  is  (2.11)  with  the  black  grid-points  at  iteration  n+1  and 
the  red  grid  points  at  iteration  n+h- 

The  SSOR-SI  and  the  SSOR-CG  method  use,  for  its  basic  iterative 
equation  the  SSOR  scheme  with  relaxation  factor  o)  to  accelerate  the 
rate  of  convergence.  The  SSOR  procedure  involves  a forward  and 
backward  sweep  of  all  grid  points  with  the  natural  ordering.  A 
symmetric  positive  definite  iteration  matrix,  «S  , is  obtained  from 
this  to-and-fro  sweep.  The  natural  ordering  is  used  since  the  opti- 
mum relaxation  factor  for  SSOR  with  the  red-black  ordering  iso.  =1, 
i.e.,  there  is  no  advantage  in  using  the  SSOR  procedure  with  the  red- 
black  ordering.  The  basic  iterative  equations  for  methods  based  on 
the  SSOR  method  is 


= rv  .(n)  4.  V . (n)  , . „(n+Ji)  „ ,(n+is) 

ij^i+l,j  ij^i,j+l  i-l,j  i-l,j  + i,j-l  i,j-l 


+ R.  . )/C  . . + (l-(a)u: 
13  iJ  13 


13  13  i+l.J  13  1.3+1  1-1.3  1-1.3  1.3-1  1.3-1 


Thin  cnn  be  written  in  matrix  form  an 
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+ k<“> 

0)  0) 


or 

(2.12) 


(n+1) 

u 


cS  u 


(n) 


k 

u) 


whore 


= <u 


iii 


k»>+k«) 


k = <i/ 

0}  0)  u u 


The  six  iterative  methods  investigated  in  this  study  apply  either 
Chebyshev  Acceleration  (Semi-iteration)  or  Conjugate  Gradient 
Acceleration  to  a basic  method  of  the  form 

(n+1)  (n)  , , 

u * y u +k 


where 


1 

k 

Jacobi  ; 

B 

, c 

Reduced  System  : 

F F 

B R 

• Vr-*- 

‘^B 

SSOR  : 

S = HI  £ 

U U)  01 

<U  k^^^ 

* Ol'^’O) 

+ k(^> 
01 

Both  the  Chebyshev  and  the  Conjugate  Gradient  Acceleration 
procedures  for  basic  methods  of  this  form  can  be  written  as 


(n+1) 


n+1 


(Y  + u^"^)  + (1  - 

tl 


f'n+l^^ 


(n-1) 


wher  e 

^(n)  _ ^^(n)  ^ ^ _ ^(n)  ^ 

Here  P^+l  ^n  acceleration  parameters  which  are  determined 

automatically  in  the  algorithms.  As  a reference  for  these  methods 
and  the  acceleration  algorithms  consult  Hageman  and  Young  [1]. 
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Detailed  descriptions  of  the  following  six  adaptive  algorithms 
are  given  in  Appendix  1. 

I.  Jacobi  Semi- i ter at  ion  (J-SI) 

II.  Compressed  Jacobi  Conjugate  Gradient  (CJ-CG) 

III.  Reduced  System  Semi- i t e r a t io n (RS-SI) 

IV.  Reduced  System  Conjugate  Gradient  (RS-CG) 

V.  Symmetric  Successive  Overrelaxation  Semi- 
Iteration  (SSOR-SI) 

VI.  Symmetric  Successive  Or er r e 1 axa t ion  Conjugate 
Gradient  (SSOR-CG) 

We  should  note  that  procedures  based  on  the  SSOR  method  require 
twice  as  much  work  per  iteration.  Also,  that  the  J-CG  method  re- 
quires exactly  twice  the  number  of  iterations  as  does  the  RS-CG 
method . 
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3 . ITPACK  Structure  and  Use 

The  ITPACK  collection  of  codes  preforms  various  tasks  which 
are  accomplished  in  individual  modules.  The  basic  modules  are  (1) 
grid  definition,  (2)  generation  of  the  nonzero  coefficients  of 
the  linear  system,  (3)  definition  of  the  ordering  vector  for  the 
grid  points,  (4)  initialization  of  the  unknown  vector,  (5)  solu- 
tion by  an  iterative  method,  and  (6)  output  of  results. 

The  grid  definition  is  accomplished  by  the  subroutine  REGION. 

In  its  present  state,  REGION  accepts  a polygonal  parameterization 
of  the  domain  of  interest.  However,  this  parameterization  must 
be  established  using  horizontal,  vertical,  and  forty-five  degree 
lines.  Consequently,  it  is  only  designed  to  accept  uniform  mesh  spacing. 
It  will  however  allow  regions  with  holes  in  them.  REGION  generates 
a rectangular  grid  and  defines  an  integer  array  GTYPE  such  that  for 
each  grid  point  (I,J)  the  value  of  GTYPE  is  either  1,2,  or  3 which 
indicates  either  interior,  boundary,  or  exterior  grid  points, 
respectively.  REGION  also  defines  arrays  GRIDX  and  GRIDY  which 
contain  the  coordinates  of  the  grid  points  in  the  x and  y direction. 

In  addition,  REGION  defines  the  minimum  and  maximum  x and  y values 
(AX , BX , AY , BY) , the  actual  number  of  grid  points  in  each  direction 
(NGRIDX, NGRIDY) , and  the  total  number  of  grid  points  (NGRPTS) . The 
remainder  of  the  ITPACK  code  needs  only  the  grid  information  gen- 
erated by  REGION  and  not  the  parameterization.  A complete  listing 
of  REGION  is  given  in  Appendix  3 along  with  additional  details 
on  the  use  of  this  subroutine. 

The  next  task  is  that  of  generating  the  nonzero  coefficients 
of  the  associated  linear  system.  This  is  accomplished  in  the 
Fortran  module  FIVEPT  which  is  currently  designed  to  handle  only 
self-adjoint  elliptic  operators.  Therefore,  the  linear  system  is 
symmetric  and  a symmetric  storage  scheme  can  be  used.  These  non- 
zero coefficients  are  placed  in  a four-column  array  COEF  as  follows: 


C0EF(IJ,1) 
C0EF(IJ,2) 
C0EF(IJ,3) 
C0EF(IJM1, 2) 


center  coefficient  at  (I,J) 
north  coefficient  at  (I,J) 
east  coefficient  at  (I,J) 
south  coefficient  at  (I,J) 
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C0EF(IM1J,3)  = west  coefficient  at  (I,J) 

C0EF(IJ,4)  « right-hand  side  at  (I,J) 

where 

IJ  = I + (J-1)*NGRIDX 

IJMl  = I + (J-2)*NGRIDX 

IMIJ  = (I-l)  + (J-1)*NGRIDX 

The  basic  iterative  equation  then  becomes 

U(IJ)  » (C0EF(IJ,3)*U(IP1J)  + C0EF(IJ,2)*U(IJP1) 

+ C0EF(IM1J,3)*U(IM1J)  + C0EF(IJM1,2)*U(IJM1) 

+ C0EF(IJ,4))/C0EF(IJ,1) 

where 

IPIJ  = (I+l)  + (J-1)*NGRIDX 
IJPl  = I + J*NGRDIX  . 

FIVEPT  requires  the  subroutine  PDE  which  is  user  supplied  or 
generated  by  the  ELLPACK  control  program.  PDE  computes  the  coeffi- 
cients of  the  self-adjoint  elliptic  operator  at  the  point  (x,y). 

A sample  of  the  use  of  PDE  is  given  in  Appendix  2. 

To  allow  extensions  to  three-dimensional  problems,  a one- 
dimensional array  is  used  for  the  unknown  vector  with  the  elements 
ordered  so  that  a linear  sweep  through  this  array  is  the  same  as 
proceeding  through  the  grid  points  with  the  natural  ordering.  At 
present,  four  orderings  have  been  coded  and  tested,  namely,  the 
natural  ordering  (NATORD) , the  red-black  ordering  (RBORD) , the 
diagonal  ordering  (DIAGORD) , and  the  spiral  ordering  (SPIRORD) . 

Each  of  these  subroutines  defines  the  arrays  NDXEQ  and  INVNDX. 

NDXEQ  is  defined  in  such  a way  that  J = NDXEQ(I)  means  that  the  I-th 
point  that  is  swept  is  actually  the  J-th  point  in  the  natural 
ordering.  INVNDX  is  the  inverse  index  array  defined  such  that 
INVNDX (NDXEQ ( I ) ) ■ I.  This  convention  is  outlined  in  [ 5 ] . These 
arrays  enable  the  same  code  to  use  any  ordering  specified.  It  is 
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i.r.  t c r e s D i riii  to  note  that  for  efficiency  in  production  software  one 
might  want  to  write  code  so  that  the  ordering  was  encoded  into  the 
iterative  algorithm;  however,  this  would  not  allow  any  versatility. 

The  next  ITPACK  task  is  to  initialize  the  unknown  vector  in 
the  subroutine  INTUNK  which  uses  the  user  supplied  (or  ELLPACK 
generated)  routines  APXUNK  and  BCOND.  The  subroutine  APXUNK  com- 
putes the  initial  appr oxima tion  (or  guess)  for  the  unknown  (or 
solution)  vector.  When  no  information  is  available  the  value  of 
zero  is  taken  for  the  initial  guess.  The  subroutine  BCOND  computes 
the  values  of  the  boundary  grid  points.  Subroutine  INTUNK  sets 
the  elements  of  the  array  UNKNWN,  which  corresponds  to  interior 
grid  points,  to  the  values  supplied  by  subroutine  APXUNK.  Other 
elements  which  correspond  to  boundary  grid  points  are  set  to  values 
supplied  by  BCOND  while  exterior  grid  points  are  set  to  zero. 

Certain  input  data  besides  the  parameterization  information 
for  REGION  and  the  subroutines  PDE  APXUNK  and  BCOND  must  be 
supplied.  The  following  is  a list  of  all  the  necessary  input  data; 


LEVEL 

NGRDXD 

NGRDYD 

MXNCOE 

MXNEQ 

ITMAX 


ZETA 

EPSI 


Controls  output  of  REGION 

Maximum  number  of  grid  points  in  the  x direction. 
Also  the  first  dimension  of  GTYPE,  and  dimension 
of  GRIDX. 

Maximum  number  of  grid  points  in  the  y direction. 
Also  the  second  dimension  of  GTYPE  and  dimension 
of  GRIDY. 

Set  equal  to  4 for  five-point  stencil.  Later 
a value  of  6 will  indicate  a nine-point  stencil. 
Dimension  of  UNKNWN  and  first  dimension  of  COEF. 
Commonly  taken  to  be  NGRDXD*NGRDYD . 

Upper  bound  on  number  of  iterations  the  user 
will  allow  the  method  to  take  before  convergence. 
If  ITMAX  is  reached,  the  methc  i will  stop  and 
exit  naturally.  Note  that  the  stopping  criteria 
may  not  be  satisfied. 

Tolerance  level  in  stopping  test  (usually  10  ^) . 
Tolerance  level  in  root  solving  and  checks  in 
division  by  zero  . (usually  10  ^). 
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CASE 


Initial  guess  of  largest  eigenvalue  of  the 
iteration  matrix.  If  no  information  is 
known,  CME  = 0.0  is  acceptable. 

Initial  guess  of  smallest  eigenvalue  of  the 

iteration  matrix.  If  no  Information  is  known 

{0.0  if  CASE  = false 

-1.0  if  case  = true 

A logical  variable  to  indicate  which  case  of 
the  adaptive  procedure  is  used. 

A factor  used  in  the  adaptive  procedure 
(usually  F = .75). 


A workspace  area  must  be  supplied  in  blank  common.  The  size 
of  this  workspace  varies  for  each  method  and  the  variable  MXNEQ. 

The  workspace  is  used  in  various  capacities,  but  is  primarily 
needed  for  the  auxiliary  storage  utilized  in  the  Iterative  algorithms 
At  present,  the  workspace  array  WORKSP  must  be  dimensioned  as 
follows  for  each  iterative  method: 


Minimum  Value  of  Dimension  for 


Method 

J-SI 

CJ-CG 

RS-SI 

RS-CG 

SSOR-SI 

SSOR-CG 


WORKSP 
3*MXNEQ 
3*MXNEQ  + 200 
2*MXNEQ 
4*MNXEQ  + 200 
5*MXNEQ 
6*MXNEQ  + 200 
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4.  Test  Problems  and  Regions 

In  order  to  test  the  code  written  to  date  for  ITPACK,  six 
test  partial  differential  equations  with  known  solutions  and  six 
reliant,  were  selected.  The  test  cases  were  desif^ned  so  that  the 
behavior  of  the  six  iterative  algorithms  could  be  monitoreu. 

The  test  equations  cover  a wide  range  of  self-adjoint  operauor 
of  the  form 

L (u)  = f 

For  each  of  the  test  problems,  u is  known  over  a region  R. 

Furthermore,  on  the  boundary  of  R,  the  function  u is  set  to  the 

true  solution  of  the  problem.  For  each  iterative  method,  the 

initial  approximation  for  u on  the  iterior  of  R was  selected 

to  be  identially  zero.  , - 

2 3 ^ 3 ^ 

In  the  following  equations,  V • » j • + r • and 

3 3 

y.  > -2-  . + _ . 

3x  3y  • 

The  test  problems  are  as  follows: 

(1)  V^u  *»  f 


where 


f = bxye^"*"^  (xy+x+y-3  ) 
“ 3xye^^^(x-l) (y-1) 


(2)  (e  + (e  ^^y)y  “ u/(l+x+y)  = f 


where 


= tt{x  sinCtrx)  cos(7ry)  + 3ye^’'^  cos(t;x)  sin(iTy)} 

+ sin(Trx)  sin(iry)  { (2y^-Tr^)  e^^^  - it^  - e*^ / ( 1+x+y ) } 


true 


e^^sin(Try)  sinCfrx). 


(3) 


where 


V.[  (1  + sin(  ^ (x+y)  Vu]  = f 


f = 8 [1  + sin ( ^ (x+y) ) ] 


true 


+ TT^sin(  j (x+y)  ) [ (x  - ^)  ^ + (y  - 
2[(x--|-)^  + (y--|)^]/[l  + sin(^  (x+y ) ) ] 


(4) 


where 


V^u  = f 


2 2 

f = 8(x  +y  -x-y) 


Ufrue  ° 4xy (x-1) (y-1) 


(5) 


where 


V u - lOOu  = f 


true 


300  cosh ( 20y) / cosh ( 20 ) 

cosh ( lOx) /cosh ( 10)  + cosh ( 2 Oy ) / c o sh ( 20 ) 


(6)  (A(x)u^)^  + (C(y)Uy)  = f 


where 
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(2+x)  e^  - ( 1+y  ) sin  (Try ) + cos(TTy),  x,y£(0,is] 
( 1-x)  e’^  - TT^  ( 2-y  ) sin  (Try ) - cos(TTy),  x,ye(4tll 


u =e  +sln(TTy) 

true 


and  where 


For  each  of  the  six  test  regions,  the  vertices  indicated  by 
dots  are  as  follows: 

(1)  (0,0)  (1,0)  (1,1)  (0,1) 

(2)  (0,0)  (1,0)  (1,1)  (.8,1)  (.8, .8)  (.6, .6)  (.4, .6)  (.2, .8) 

(.2,1)  (0,1) 

(3)  (0,0)  (1,0)  (1,.4)  (.6, .4)  (.4, .6)  (.4,1)  (0,1) 

(4)  contour  1 - (0,0)  (1,0)  (.5, .5)  (.5,1)  (0,1) 

contour  2 - (.2, .2)  (.2, .4)  (.4, .4)  (.4, .2) 

(5)  (.2,0)  (.6,0)  (1,.4)  (1,.6)  (.6,1)  (.6,1.2)  (.2,1.2)  (0,1) 

(0,.8)  (.2, .6)  (0,.4)  (0,.2) 

(6)  contour  1 - (.6,0)  (1,.4)  (1,.7)  (.7, .7)  (.7, .9)  (.4, .9)  (.4, .6) 
(0, .6) 

contour  2 - (.  4 , . 3)  (.  4 , . 5)  (.  45, . 5)  (.45  . .35)  (.  5 , . 35) 

(.  5 , . 5)  (.  55  , . 5)  (.  55  , . 3) 

contour  3 - (.  7 , . 3)  (.  7 , . 45)  (.  65, . 45)  (.  65, . 5)  (.  8 , . 5)  (.  8 , .45) 
(.  75  , . 45)  (.  75, . 3) 


The  Identification  grid  array  GTYPE  generated  by  subroutine 
REGION  is  given  in  Appendix  3 for  each  of  these  regions. 
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1 

5 . Numerical  Results 

In  this  section  we  will  discuss  the  results  of  numerical 
test  runs  with  the  ITPACK  code.  The  first  set  of  test  runs 
were  on  the  unit  square  with  h = 1/40.  This  was  done  to  compare 
the  results  of  the  six  iterative  methods  on  a common  region  over 
a variety  of  problems.  All  six  methods  were  run  on  the  six  test 
equations  described  in  Section  4 with  the  following  initial  data:  \ 


F = 

.75 

CME  = 0.0, 

EPSI  = 

.000001, 

SHE  = 0.0, 

ZETA  = 

. 000001, 

CASE  = .FALSE. 

The  red/black  ordering  was  used  with  J-SI,  RS-SI,  RS-CG,  and  CJ-CG. 
The  natural  ordering  was  used  with  SSOR-SI  and  SSOR-CG.  All  runs 
were  made  on  a CDC  6600  with  the  MNF  compiler  and  UT2D  operating 
system. 

For  the  second  set  of  test  cases  we  considered  the  elliptic 
operator  described  as  test  equation  (2)  in  Section  4.  Each  itera- 
tive method  was  used  with  the  five  test  regions  described  in 
Section  4 with  h=l/20.  All  input  data,  initial  conditions,  etc. 
used  were  the  same  as  those  selected  in  the  first  set  of  test  cases. 

The  following  tables  represent  these  runs  and  give  the  resulting 
data  for  comparison.  Each  block  of  the  tables  has  the  form 


where 


20 


A * Number  of  iterations  until  convergence 
^ ^^'^true  ^computed^^  ^ II  ^true^^ 


C 


For  J-SI,  RS-SI,  SSOR-SI,  SSOR-CG,  a list  of 

iterations  numbers  where  parameters  were  changed 
For  CJ-CG  and  RS-CG,  the  iteration  number  where 

a new  estimate  of  CM£  was  last  calculated 

V. 


D “ Last  estimate  of  CME. 

For  further  comparison  of  the  iterative  methods,  contour 
plots  of  error  distributions  between  computed  and  true  solutions 
were  generated.  All  of  these  test  cases  used  the  self-adjoint 
operator  described  as  problem  (4)  on  the  unit  square  with  h“l/2C. 
Problem  (4)  was  used  so  there  would  be  no  discretization  error 
from  the  five  point  difference  equations.  The  contour  plots  are 
of  the  function  z(x,y)  defined  as 


|u^  (x,y)  -u  ^ j(x,y)| 

z(x  y)  •»  true  computed  ' 

SCALE 

where 

SCALE  - max  I “true  " “computed  I ' 

Figures  1 thru  4 are  the  error  distributions  at  convergences 
of  SSOR-SI,  SSOR-CG,  RS-SI,  and  RS-CG  respectively.  Figures  5 and 
6 show  the  error  distribution  of  RS-CG  after  five  and  ten  itera- 


tions, respectively.  The  scaling  factor,  SCALE,  which  is  the  max- 
imum pointwise  absolute  error  is  given  for  each  case. 


Table  1.  Test  Problems  Over  the  Unit  Square  With  h-1/40 


Table  2»  Test  Problem  (2)  With  h = 1/20 


SSOR-CG  method  error  distribution  at  convcryoncc 


li'.ur  o 


(SCALE  = 1.583368  x 10 
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Description  of  Adaptive  Procedures 
(Equations,  Flow  Chart,  and  Algorithm) 

Precise  descriptions  of  the  following  six  adaptive  iterative 
algorithms  are  stated  in  this  appendix. 

I.  Jacobi  Semi-iteration  (J-SI) 

II.  Compressed  Jacobi  Conjugate  Gradient  (CJ-CG) 

III.  Reduced  System  Semi-iteratl  on  (RS-SI) 

IV.  Reduced  System  Conjugate  Gradient  (RS-CG) 

V.  Symmetric  Successive  Overrelaxat ion  Semi- i t era 1 1 on 
(SSOR-SI) 

VI.  Symmetric  Successive  Overrelaxation  Conjugate  Gradient 
(SSOR-CG) 


For  each  method,  a list  of  equations,  a flow  chart,  and  an  algorithmic 
description  is  given.  The  latter  description  details  exactly  the 
adaptive  procedure  used  in  the  ITPACK  code.  The  mathematical  deriva- 
tion for  each  of  these  methods  can  be  found  in  Hageman  and  Young  [Ij. 
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I . J-SI;  Jacobi  Seml-lteratlve  Equations 


(1)  Adaptive  Parameters 

V . 2/  (2-m^-,.^)  , / (2-M^-„^) , 

t . (I  - [l-cr|l*’)/{l  + (l-(rj)^) 


(2)  Acceleration  Parameters 


n+1 


l/[l-a^/2],  n = s+1 

hi 


l/fl-(o-„/2)  p ],  n > s+1 
t n 


(3)  Residual  Vector 


.(n) 


+ c - 


n = s+1 

n > s + 1 


(4)  Iteration  Vector 


(s+1)  _(s)  (s) 

u =y5  +u  , 


n = s 


f »(“)  (v)  \ /I  \ (n-1) 

‘ = Pn+1^Y&^  + M + (1-p^^pu^  ^ 


(5)  Stopping  Test 
d(n)  = 6^*^^ 

T 

STEST  = [1/(1-M^) ] [d(n)/(u^”^ 
If  STEST  < then  exit. 


n > s+1 


J-SI  (Continued) 
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6.  Changing  Parameter  Test 

QA  = [d(n)/d(s) 

QT  = 2r^"‘^^'^^/(l+r^"'®^) 

F 

If  QA  > QT  , then  change  parameters. 

7.  Rayleigh  Quotient  Vector 

~(n+l)  ^ ^(n)  ^ g(n) 

+ c - 

8.  Computing  new  M^.  and  m^, 

Z = (l+r^""®^{QA  + [QA^-QT^1^^^)/2 
X = 

<J  = (X+r/X)/(l+r) 

if  n = 0 


Mi  = 


J-v 

[Me^>  0(2-11^-10^)  ] /2,  otherwise 

(~(n+l)  ^ ^ 

Mg  = max{Mj^,M2} 


«2  = 


if  case  I 
if  case  II 


r 


not  changed,  if  case  I 


, if  case  II 


3A 


J-SI;  Jacobi  Seml-lterat Ive  Algorithm 


Input  ,CME,SME,F,EPSI,ZETA,CASE,ITMAX,LEVEL 

Set  N:=0,S:=0 

T - 1 

Compute  CNRM:=c  *D*c,  where  c=D  b 

If  CNRM  < EPSI,  then  go  to 


EXIT 


I START  I If  N > ITMAX,  then  go  to  (eXIT 

Tf  M^CJ.1  l-Kn« ^ — D ^ ) _l_„ 


If  Ni^S+I,  then  compute  6 =B*u  ^^^+c-u 

else  compute  6 ^ ^ =GAMMA  • 6 ^ ^ + (1-GAMMA)  • <S  ^ ^ 


(Test  for  stopping)  ^ 

Compute  UNRM:=u^"^  *D*u^”^ 


DELNRM: =6 


*D*  6 


(n) 


If  UNRM<CNRM/2,  then  set  UNRM:=CNRM/2 

L 

Compute  STEST:=(DELNRM/UNRM)^/ (1-CME) 
If  STEST  < ZETA,  then  go  to  IeXItI 

(Test  for  changing  parameters) 


Tf  N = 0,  then  go  to  IcHANGeI 

u 

Compute  QA:=(DELNRM/DELSRM)^,  P:=N-S, 

P / 2 P 

QT:=2R  ' /(1+R  ) 


If  QA>QT  , then  go  to  |CHANGE'| 


(Preform  iteration  with  current  parameters) 
If  N=S+1,  then  compute  RHO : = 1 . 0 / ( 1- SIGE^ / 2 ) 
else  compute  RHO : =1 . 0/ ( 1-RHO • S IGE^/ 4 ) 
Compute  Cl :=RH0. GAMMA, C2:RHO,C3=l-RHO, 
^(n+l)^Ci. 6(")+c2'u^"^+C3-u^"“^^ 


Go  to  lENDITl 


r 


J-SI  (continued) 
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IcHANGeI  (Change 
Compute 


parameters) 

~ (n+l)^^(n)^g(n) 


If  N=0,  then  set  ZM1:=CME, 

else  compute  Z : = (1+R^) (QA+(QA^-QT^)^) / 2 , 

X:=zl/P 

SIG: = (X+R/X) / (1+R) 
ZM1:=CME+SME+SIG- ( 2-CME- SME ) ) / 2 


If  CASE=.TRUE.  , then 
else  compute  ZM2: 


compute  ZM2:=6^"^  *0*6^"''’^^ 

= *D*6  VdELNRM 


/DELNRM 


Set  CME:=max{ZMl,ZM2} 

If  CASE=. FALSE. , then  set  SME;=-CME 
Compute  S1GE:= (CME-SME) / (2-CME- SME) 

GAMMA: =2/ (2- CME-SME) 

R:=(l- (1-SIGE^)^) / (1+(1+SIGE^)^) 
Set  S:=N 

DELSRM:=DELNRM 

RH0:=1 

Print  N,ZM1,ZM2,CME 
Compute  u^"’*’^^=GAMMA- 


I ENDIT I (End  of  iteration  step) 

Print  N , UNRM"^  , STEST  , QA,  QT^  , CME  , RHO  , GAMMA 
Set  N:=N+1 
Go  to  ISTARTl 


EXIT 


(Exit  iteration  algorithm) 

Compute  UNRM:=u^"^  *D*u^"^ 
Print  N,UNRM^ 

If  LEVEL  > 2,  print  u^"^ 
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II*  CJ-CG;  Compressed  Conjugate  Gradient  Equations 


(1)  Residual  Vector  (non- recursive  computation) 


(n)  (n) 

R R B R 

c.(n)  „ (n)  (n) 

o'  = F„u'  ' + c - u' 

B B R B B 


(2)  Acceleration  Parameters 


R B ’ 


n = 0 


n - 0 


i/[i  - (l/Pn)(dB(n)/dR(n-l))]  , n>0 


(3)  Iteration  Vector 


. 


n = 0 


+ (l-a)u^"'2\  n >0 


(4)  Residual  Vector  (recursive  computation) 


^n.2^*»R 


(n4l) 


CJ-CG  (continued) 
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(5) 


Stopping  Test 

(a)  Compute  which  is  the  largest  eigenvalue  of  the  symmetric  nxn  tri- 
diagonal  matrix  (1  < i < n) 


X /2 

(b)  STEST  = [/27  (1-M^>]  [dg(n)/(Ug(n)^*Dg*u^"^] 

If  STEST  < 8,  then  exit. 

Note;  If  < e,  then  an  acceptable  estimate 

of  M 

E 


has  been  obtained  and  computation  5a  is  omitted. 
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Compute  Residual 
Vector 

I (Non- recursively) 


Flow  Chart  2;  CJ-CG  Method 


^ Input  EPSI.ZETA, ITMAX, LEVEL 

Set  N;=0 

T 

; Compute  CNRM : =kg*Dg*kg 

; If  CNRM<EPSI,  then  go  to |eXIT 

START  If  N > ITMAX,  then  go  to  [exIT 
If  N < 4,  then  go  to  ONE 

If  [CME  - CMOLDj /CME  < EPSI,  then  go  to [TWO 

|one|  (Determine  new  CME) 

Set  CMOLD:=CME 

I If  N=0,  then  set  CME:=0 

1 

i Else  set  CME:=maxlmum  eigenvalue  of  the  trl-dlagonal  matrix 

{ [(RH0^-1)/(RH0^RH0^^^)]*^,0,  [ (RH0^_^^-1)  / (RHO^_^^RHO  . ) ] ‘'^  } , 

for  1 ^ 1 i N 


CJ-CG  (continued) 


EXIT 


If  N=0,  then  set  RH0HAT:=1, 

else  compute  RHOHAT : =1+RH0„^. (1-RH0„^, ) • ( 1-RH0„) /RHO 

N+/  N+1  N 

Compute  GAMMA: =RHO„^- • RHO„^, /RHOHAT 

N+z  N+1 

Set  Cl: =RHOHAT- GAMMA, C2 : =RHOHAT , C3 : =1 -RHOHAT 

_ ^ (N+2)  (N),^-  (N-2) 

Compute  u_  =C 1 • 6 +C 2 • u +C3*u_ 

D D O D 


Compute  V 


B 


F 

B R 


Set  c1:=RH0jj^2»  ^ ^ : =1-RH0j^^2 
Compute  6g  =C1.^^+C2•6g 

Print  N.UNRM^, STEST,CME, RHO, GAMMA 
Set  N=N+2 
Go  t o 


START 


Compute  u^*^^  =F*Ug^^+C^ 


UNRM: =u 
Print  N,UNRM'^ 


(N) 


*D*u 


(N) 


IF  LEVEL  > 2,  then  print  u 


(N) 


1 


III.  RS-SI;  Reduced  System  Semi-Iterative  Equations 


(1)  Adaptive  Parameters 

Y = 2/(2-m|), 


''e  = 


r = {1  - [1-M^]^)/{1+ 


(2)  Acceleration  Parameters 


1 

n+1 


(3)  Residual  Vector 

r 


l/[l-a^/2] 


1/[1- (0-^/2)  ^p^]. 


n = s+1 


n > s+1 


(n)  _ (n) 

= F„u;  ' + c„ 
R R B R 


-(n)  (n)  (n) 

o =Fu  +c  - u 
B B R B B 


(4)  Iteration  Vector 


(s+1)  As')  (s) 


» P.,,(y6 


(n)  (n)  , . (n-1) 

wr-B  +"b  ' * - 


(5)  Stopping  Test 


B^  B B B 


STEST  = [ n/2  /(1-M^  ] [dg(n)  /(Ug(n)^*Dg*Ug(n)  ] 


1/2 


If  STEST  < ^ then  exit. 
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> s + 1 


1 
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RS-SI:  Reduced  System  Semi-Iteration  Algorithms 


f I START 


Input  u^®\cME,F,EPSI,ZETA,  ITMAX 
Set  N:=0,S:=0 
Compute  CNRM:=k3*D  *k 

ODD 

If  CNRM<EPSI,  then  go  to  (eXIT  | 


IF  N > ITMAX,  go  to  (exit! 

(N)  (N) 

Compute  u^  " ^R*^B  ‘^R 


(N) 


+ c, 


u 


(N) 

B 


(Test  for  stopping)  ^ 

Compute  UNRM:=uf^^  *D„*u^^^ 

B ^ B B 

DELNRM:=6^^^ 

B B B 


RS-SI  (continued) 


I CHANGE  i 


(Change  parameters) 

If  N«0,  then  set  ZM1:=CME, 

else  compute  Z : • ( 1+R^^ ) (QA+ (QA^-QT^ ) *^)  / 2 

X:.zl/(2P) 


ZMl :=(X+R/X)/ (1+R) 


V 


R 


: = F*6 
R 


(N) 

B 


ZM2:=  (v^*DJ^*VJ^^ELNRM)^ 


Compute  CME  : =*max{  ZMl , ZM2  } 


2 2 

Compute  SIGE:=CME  /(2-CME  ) 

GAMMA:=2/ (2-CME^) 

R:  = (1-  (l-CME^)*^)  / (1+(1-CME^)*^) 


Set  S : =N 

DELSRM: =DELNRM 
RH0:=1 

Print  N,ZM1,ZM2,CME 
Compute  u ^ =GAMMA • 


6<^Uu 


(N) 


ENDIT 


Print  N.UNRM^, STE ST, QA.QT^.CME.RHO, GAMMA 
Set  N:=N+1 
GO  to  ISTART I 


(N)  (N) 

Compute  uj  ^=F*u; 

(N)T 

UNRM:=u^  ’ *Dj^*u 
PRINT  N.UNRM** 


If  LEVEL  > 2 , print  u 


(N) 


(N) 

R 
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IV.  RS-CG;  Reduced  System  Conjugate  Gradient  Equations 


(1)  Residual  Vector  (non- recursive  computation) 


u^“^  = F *u^"^  + c 
R R B R 


5^^)  = F . c - u^"> 

B B R B B 


(2)  Acceleration  Parameters 

,(n) 


V = F *5' 
R R B 


V = F *v 
B B R 


d (n)  = 

B B B B 


Vl  = .D^«  )/d^(„)] 


n =r  0 


n+l 


(3)  Iteration  Vector 


(n+1)  ^ 
“b 


„(n)  (n) 

" “b 


n = 0 


Pn.l  ^ ^ « >0 


.(n-1) 
"n+l'^'B 


(4)  Residual  Vector  (recursive  computation) 

f ^rB " 


,(n+l) 


, n = 0 


Ip  1 (v  i'^o  ' ' (^"P 

•^n+l  nfl  B n^l  B n^  1 


3 (n-1) 
B 


n > 0 


J 
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RS-CG  (continued) 

(5)  Stopping  Test 


(a)  Compute  which  is  the  square  root  of  the  largest  eigenvalue 

of  the  symmetric  nxn  tridiagonal  matrix  (s+1  < i < n) 


) U.L.] 

LVviPj.iVjPi.;’  I ''i / ’ I Vi ViVi ' 


(b)  STEST  = [ /2/(l-M5][d  (n)/(u^"^  *D 

E S B B B 


If  STEST  < p,  then  exit. 


Note;  If  |m^  ^ then  an  acceptable  estimate 

of  Mg  is  available  and  computation  5a  is  omitted. 
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RS-CG;  Reduced  System  Conjugate  Gradient  Algorithm 


Input  u'  \EPSI,ZETA, UMAX, LEVEL 
Set  N:=0 

Compute  CNRM: =k3*D*k„ 

If  CNRM<EPSI,  then  go  to  EXIT [ 


If  N > UMAX,  then  go  to  | EXIT| 

If  N < 4,  then  go  t o | ONE | 

If  I CME-CMOLdI /CME  < EPSI,  then  go  to  fTWO 


(Determine  new  CME) 

Set  CMOLD:=CME 

If  N=0,  then  set  CME:=0 

else  set  CME=square  root  of  the  maximum  eigenvalue  of  the 
tridiagonal  matrix 

{[ (RHO  -1) / (GAMMA. • RHO  . -GAMMA. -RHO)  (1-1/GAMMA.) , 

i Ill  1 

[ (RHO^^^-l)/(GAMMA^-RHO^-GAMMA^^^ -RHO^^j) 

(Test  for  stopping) 

(N)^  (N) 

Compute  UNRM:=u„  *D*uy 
B B 

DELNRM:=6^^^ 

B B 

If  UNRM < CNRM,  then  set  UNRM:=CNRM 

Compute  STEST:  = (2-DELNRM/UNRM)'^/  (1-CME^) 

If  STEST  < ZETA,  then  go  to  EXIT 

(N) 

Compute:  " ^B*^R*^B 


GAMMAjj_^j^:=l/(l-6^‘ 


*D„*v_A)ELNRM) 
B B 


If  N=0,  then  RHO^:=0 

else  compute  RHO„^, : =1 / ( 1-GAMMA„^, DELNRM/ (GAMMA„DELSRM • RHO  ) ) 
^ N+1  N+1  N N 


Set  DELSRM:*DELNNM,Cl:-RHO^^^-GAMMAj^_^^^,C2:«RHOj^^j^,C3:*l-RHO^^^, 


C4  : ( 1-GAMMA^, ) 

N+1  N+i 


! 


RS-CG  (continued) 
(N+1) 


Compute:  u 


B 


Cl*6g^^+C2*Ug**^+C3*Ug*^"® 


Print  N,UNRM^,STEST,CME,RHO„^, ,GAMMA„^, 

N+i  N+i 

Set  N=N+1 


Go  to  START 


I EXIT  I Compute  u:"^F*u:‘' '+C 


R 


(N) 

B ‘ *R 
T 


UNRM:=u^^^  *D*u^^^ 


Print  N.UNRM^ 

If  LEVEL >2,  then  print  u 


(n) 


I END 
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V.  SSOR-SI;  Symmetric  Successive  Over relaxation  Semi-iterative  Equations 


(1)  Adaptive  Parameters 


SPECR  = S(.8  ),  Y = 2/(2-SPECR),  a_  = SPECR/(2-SPECR) 

Cl'  b 


r = {1  - [1-SPECR^]^^^)/{1  + [1-SPECR^]^^^) 


(2)  Acceleration  Parameters 


1/ [1  - a /2]  , n = s+1 

E 


l/[l-(a^/2)  p ],  n > s+1 
E n 


(3)  Difference  Vectors  and  Residual  Vector 


V = a! 

CL)  CD 


5^"^  = <M  V + 


ui  ai 


(4)  Iteration  Vector 


(n+1)  r .c.(n)  (n)  i s (n-1) 

u = M + (l-p„^pu' 


(5)  Stopping  Test 


d(n)  = *D*A^"^ 


T 

STEST  = [(2-<jL')/(tu(l-Mg)  (1- SPECR)  ^ ]^^^[d(n) /(u^^^  *D*u^*^^  ) ] 


If  STEST  < ^ then  exit. 


SSOR-SI  (continued) 
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(6)  Changing  Parameter  Test 


QA  = [d(n)/d(s)] 


1/2 

/(l+r^“-“^) 


If  QA  > QT  , then  change  parameters. 


(7)  Computing  new  S'  and  M 


Z = (QA  r [QA^-QT^]^^^}/2 

l/(n-s) 


X = Z 
0-  = (X  + r/X)  /(1+r) 


SPECR 


, if  n = 0 


[SPECR  + o-(2-SPECR)  ]/2,  otherwise 


~(n+l)  (n)  _(n) 

u'  ^ = u'  + 6 


^(n+1)  ^ ^ ~(n+l)  ^ j^(F)  _ ~(n+l) 


S = *D*i^"'^^)/d(n) 


S'  - max(SPECR,  Sj^^Sg) 


M 


[(l-S')(l-^>^  - a)(2-a)  ]/[c«(a'-l-S')  ] 

T 

V = d(n)  = 5^”^  *0*6^’^^ 

T 


(6^”^  *D*v)/d(n)  , if  case  I 


^2  = 


(v  *D*v)  /d  (n") 


1/2 


if  case  I 


= max{M^,M^,M2) 


SSOR-SI  (continued) 

(8)  Computation  of  P 


^ = max  fW  fE,  , . + N , J 
iJ  i-l,J  i-l,J 


+ S [E  . , + N J } 
iJ  i,J-l 


(9)  Computation  of  uj,  SPECR  = S(«8  ) 

a' 

If  then 

h 

CD  = 2/(1  + [1-2M^44^]^''2} 

E 

SPECR  = (2-aD-tttMg) /(2-alig) 

else 

CD  = 2/{l  + 

SPECR  = CD  “ 1 
= 2 s/p  , a = 

(10)  Computing  cd* 

CD*  = 2/{l  + [l-4p]^^^) 

(11)  Selecting  cd*  Test 

Tf  1ok(^(cd*-1))  . 

log ($ (SPECR))  - to  a*. 

(Note;  $(x)  = {l  - + ri-xj^^^).) 

(12)  Computation  of  a',  SPECR  = S(.8  ) . M 

CD  E 

O)  = CD*,  SPECR  = (D-l,  M,  = ^ p 

E 
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T 
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SSOR-SI;  Symmetric  Successive  Overr elaxat Ion  Seml-lteratlve  Algorithm 


Input  ,CME,ZETA,F, CASE, ITMAX, LEVEL, EPSI 

Set  N:=0,S:=0,OMEGCHG:=.TRUE. 


Compute  BETA:  =max{  W . . Te  . . .+N.  , ."]  + S . . Fe  . . ,+  N.  . -1} 


If  CME£4*BETA,  then  compute 

OMEGA:=2/(l+(l-2-CME+4'BETA)^) 

SPECR:= (2-2 -OMEGA+CME- OMEGA) / ( 2-CME*OMEGA) 


else  compute 

0MEGA:=2/ (1+ (1-4 -BETA) 
SPECR: =0MEGA-1 
CME:=2-BETA^ 

OMEGCHG:=. FALSE. 

Print  N,CME, OMEGA, SPECR 
Compute  CNRM:*k^*D*k 
If  CNRM<EPSI,  then  go  to  (¥xIT| 


{STARTI  If  N > ITMAX,  then  go  to  |EXIT 
(N) 


Copy  u 
Compute 


into  V. 

(F) 

V = ;£  *v  + k 

0)  u 


v=<M  *v  + k^®^ 

0)  U) 


(Test  for  stopping) 
(N)T 


(N) 


Compute  UNRM:=u'‘’'  *D*u 
If  UNRM<CNRM,  then  set  UNRM:=CNRM 
Compute  DELNRM:=A^’^^*D*A^^^ 

STEST:=  [(  (2-OMEGA) /OMEGA) (DELNRM/UNRM) / (1-CME)^/ (1- SPECR) ^ 
If  STEST  < ZETA,  then  go  to  |eXIT| 

If  N = 0,  then  go  to  | CHANGE) 

If  0MECHG-. FALSE. , then  go  to  [THREE) 

IF  BETA  i 1/4,  then  go  to  |0NE| 


SSOR-SI  (continued) 
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ONE 


CHANGE 


Compute  0MEGAS:=>2/  (1+ (1-4  •BETA)'*) 
TEMPI  :=log(«>(OMEGAS-l) ) 
TEMP2:=log(({>(SPECR)) 

where  $(X)- (1- (1-X)'*)  / (1+(1-X)'*)  . 
If  TEMP1/TEMP2  < F,  then  go  to |0N¥1 

Set  OMEGA: -OMEGAS 
SPECR:-0MEGAS-1 
OMECHG:-. FALSE. 

CME:=2-BETa'* 

S :=N 

Print  N,CME,SPECR, OMEGA 
Go  to [TWO| 

(Test  for  changing  parameters) 


Compute  QA:=(DELNRM/DELSRM)' 
QT:=2-R^^^/(1+R^) 


If  QA  > QT  , then  go  to 


CHANGE 


else  go  to [three 


P : =N-S 


(Change  parameters) 

If  N=0,  then  set  SIG1:=SPECR 

Compute  Z:  = (l+R^)  (QA+(QA^-QT^)'*)/2  i 

SIGE1:=(X+R/X)/(1+R) 

SIGl :-(SPECR+SIGEl (2-SPECR) ) /2 
T 

Compute  SIG2:=A^^^  *D*  /DELNRM 

Set  SME=max{SIGl,SIG2,SPECR} 

(Determine  new  CME , OMEGA , SPECR) 

Compute  ZMl : = ( (1-SME) (1+BETA*0MEGA^) -OMEGA (2-OMEGA) ) / (OMEGA (OMEGA- 1-SMEj 


Compute  v = B*6^”^  ^ .j, 

If  CASE-. TRUE.,  then  compute  ZM2:-(6^"^  *D*v)/(6^"^ 

T 

else  compute  ZM2  :-{^(v^*D*v)  / (6  *0*6^"^)]** 


J 
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SSOR-SI  (continued) 


Compute  CME : =max{ CME , ZMl , ZM2 } 

Set  S:=N 

If  CME  — 4BETA,  then  compute 

0MEGA:=2/  (1+(1-2CME+4BETA)*^ 

SPECR : = ( 2-2 • OMEGA+CME • OMEGA) / ( 2-CME+OMEGA) 


TWO 


else  compute 

OMEGA;=2/ (1+(1-4BETA)^) 

SPECR:=OMEGA-l 
CME:=2- BETA^ 

OMEGCHG: =. FALSE. 

Print  N, CME, OMEGA, SPECR 

Compute  R:  = (1-  (1-SPECR)^)  /(  1+ ( 1- SPECR) 

SIGE:=SPECR/ (2-SPECR) 

GAMMA:=2/ (2-SPECR) 

RH0:=1 

(Special  procedure  to  recompute  and  since  OMEGA  has 

been  changed) 


THREE 


Copy  u^"^  into  v 

(F) 

Compute  V = ' 

a' 


(n)  _ ..  ..(n) 


U) 

= v-u 

DELSRM:  = A^^^ 

V = “W,  V + k 


*D*A 
(B) 

0)  ■ "0) 

(n) 


(n) 


p (n) 

t = v-u 

= GAMMA. +u^"^ 

Go  to  (ENDIT| 

(OMEGA  has  not  been  changed) 

If  N-S+1,  then  compute  RHO : =1/ (1-S IGE^/2) 
else  compute  RHO : =1 / ( 1-RHO • S IGE^ / 4 ) 

Set  Cl: “GAMMA- RHO,c2:=RHO,C3:=l-RHO 

^ ^ (N+1)  x(N)._»  (N)_,_„_  (N-1) 

Compute  u “Cl-A  +C2'u  +C3.u 
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enditI 


EXIT 


\ 

t 

r 


SSOR-SI  (continued) 

Print  N , UNRM** , STEST  , QA , QT^  , CME , RHO  , GAMMA 
Set  N:-N+l 
Go  to [start 

(N) 

Compute  UNRM=u'  ' *0*0' 

Print  N.UNRM^ 

(N) 

If  LEVEL  > 2,  print  o'* 

end! 


1 
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VI  SSOR-CG;  Symmetric  Successive  Overrelaxation  Conjugate  Gradient  Equations 


(1)  Difference  Vector  and  Residual  Vector  (non- recursive  computation) 
w ^ £ u^"^  + 


. (n) 

A'  = V - u 


6J 

(n) 


CO  a) 


(2)  Acceleration  Parameters 

V = 5^^^  - *6^^^ 


d(n)  = A^^^  *D*A^”^ 


"'n+l  " *D*v) 


n+1 


f 


■/[■-(S') 

^ ' n n ' 


y n = s+1 

(d(n)/d(n-l))j  , n > s + 1 


(3)  Iteration  Vector 


(n+1) 


r 


Y ,5^"^  + U<"> 
'n+1 


= 


f -(n)  (n)  (n-1) 

Vl^Vl^  ^ + (1-Pn+l^^ 


n = s+1 


n > s + 1 


SSOR-CG  (continued) 
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(4)  Difference  Vector  and  Residual  Vector  (recursive  computation) 

, n = s+1 


.(n+1) 


A W 

A - Y ,V 
'n+1 


{A^”^-Y„^iv}  + n > s+1 


n+1  ^n+1 


V = 6 


(n) 


V = “W  V 

0) 


- V 

(i) 

(=  S 

0) 


^ /I  \ r(^) 

Vi'^  " ’ 


,(n+l) 


n = s + 1 


Pn+l^Vl'^  ^ " > 


(5)  Stopping  Test 

d(n)  = A^’^^  *d*A^'^^ 

1/2  . .T  („s  1/2 

STEST  = [(2-co)  /(CD(I-Mg)  (1-SPECR)  1 J [d  (n)  /(u'‘  *D*u  '•'')] 

If  STEST  < then  exit. 


(6)  Changing  Parameter  Test 

= -log[®(SPECR)/'I>(SPECR/S')  ] 

\ = -log[<t(S')] 

If  i\/\)  < F,  then  change  parameters. 


s + 1 
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(7a)  Computation  of  S' 

Compute  S'  which  is  the  largest  eigenvalue  of  the  symmetric  nxn  tri- 
diagonal matrix  (s+1  < i < n) 


(7b)  Compute  new  M^, 


= I(l-S')(l+^^  - a;(2-ai)]/fa>(a)-l-S')] 

V = d(n)  = 6^"^ 

I (6^**^  *D*v)/d(n),  if  case  I 

«2  = ^ 

' T , II/2 

I (v  *D*v)/d(n)J  , if  case  II 

Mj,  = max{Mg,M^,M2} 


(8) -(12)  Same  as  SSOR-SI 


Flow  Chart  6a:  SSOR-CG  Method 
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SSOR-CG:  Symmetric  Successive  Overrelaxation 
Conjugate  Gradient  Algorithm 


Input  u^°\cME,  F.ZETA, CASE, UMAX, LEVEI,  EPSI 
Set  N:-0, S:=0,0MEGCHG:=.TRUE. 

compute  BETA:=max{W^.[E^_^.+N^_y]+S^j[E^._^+N^._J} 

If  CMEIA’BETA,  then  compute 

OMEGA:=2/  (l+(  1-2-  CME+4- BETA) 

SPECR: = (2-2 • OMEGA+CME • OMEGA) / ( 2-CME • OMEGA) 


else  compute 

OMEGA:  =2/  (l+(  1-4  BETA)*^) 
SPECR: =OMEGA-l 

h 

CME:=2-BETA^ 


OMEGCHG:=. FALSE. 


Print  N,CME, OMEGA, SPECR 
Compute  CNRM:=k^*D*k 


If  CNRM<EPSI,  then  go  to  |EXIT; 
Copy  u^^^  into  v. 

(F) 

Compute  V = £ *v  + k 


CL- 


- v-u^°> 

V = V.  *v  + k^®^ 
U)  U) 


c(0) 


v-u 


(0) 


I startI 


If  N >ITMAX,  then  go  to  |eXIT 


(Test  for  stopping) 

T 

Compute  UNRM:=u^”^  *D*u^"^ 

If  UNRM<CNRM,  then  set  UNRM:»CNRM 
Compute  DELNRM:=A^"^*D*A^"^ 

STEST  : = [_  (2-OMEGA)  / (OMEGA)  (DELNRM/UNRM)  / (1-CME)]  (1-SPECR) 


T 


SSOR-CG  (continued) 
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If  STEST  < ZETA,  then  go  to  |EXIT| 

If  OMEGCHG=. FALSE. , then  go  to  |THREe1 
If  N = 0,  then  go  to  |THREE| 

Else  set  SME : =maximum  eigenvalue  of  the  tri-diagonal  matrix 
{[  [1-1/Y^]  , 

[ (Pi_^l-1)/ (Yi+iPi+lYiPi)  ]^}  , for  s + l<l<N 


If  BETA  >1/4,  then  go  to  |0NE| 

Compute  OMEGAS  : =2/  (l+(  1-4-  BETA) '^) 
TEMPI :=log($ (OMEGAS- 1) 
TEMP2:=log(<l>(SPECR)) 
where  $(X)  = (1- (1-X)^) / (1+(1-X)^) 


If  TEMP1/TEMP2  < F,  then  go  to |oNe| 
Set  OMEGA: =OMEGAS 

SPECR:=0MEGAS-1 
OMEGCHG:=. FALSE. 

iff 

CME:=2*BETA^ 

S : =N 

Print  N, CME, SPECR, OMEGA 
Go  to  I TWO I 


(Test  for  changing  parameters) 
If  SPECR  > SME,  then  go  to 
Compute 


THREE 


•V  _ 1 „ I ^(SPECR) 

^1  I /SPECR> 


4> 


\ SME 


X2  = -log($(SME)) 


If  then  go  to 


THREE 


SSOR-CG  (continued) 
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(Determine  new  CME , OMEGA , SPECR) 

Compute  ZM1;=((1-SME) (1+BETA • OMEGA^) -OMEGA • (2 -OMEGA)) 

/ (OMEGA(OMEGA-l-SME) ) 

Compute  „ 

r T T 

If  CASE=.TRUE.,  then  compute  ZM2=(6^"^  *D*V)/(6^”^  *D*6^"^) 

T 

else  compute  ZM2=  [|(v^*D*v)  / (6  *D*6^^^)J^ 

Compute  CME:=max{CME,ZMl,ZM2} 

Set  S:=N 

If  CME  < 4*BETA,  then  compute 

OMEGA:=2./<l+(l-2-CME+4-BETA)^) 

SPECR: = (2-2 •OMEGA+CME* OMEGA) /( 2-CME • OMEGA) 
else  compute 

OMEGA: =2. / (1 . + (1 . -4 • BETA)^) 

SPECR: =OMEGA-l 
u 

CME:=2.BETA’^ 


I TWO 


CMEGCHG:=. FALSE. 

Print  N, CME, OMEGA, SPECR 

(Special  procedure  to  recompute  and  since  OMEGA  has  been 

changed) 

Copy  u^”^  into  v 

(F) 

Compute  V = + kj; 

= v-u^"^ 

V = <«  V + 

0)  U) 

= v-u^"^ 


THREE 


DELNRM  = A^'^^  *D*A^"^ 


Copy  A into  v 
Compute  V = 

V = A-v  ^ 

'*'n  + 1 ^ DELNRM/(A^"^  *D*v) 


^ N+l“^ 


If  N=S,  then 

else  compute 

: =1/ (I-Yjj+jDELNRM/ (YfjP^DELSNM)  ) 


68 


Appendix  2 
Sample  Problem 
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The  current  ITPACK  routines  can  best  be  explained  by  looking 
at  the  code  for  a sample  problem.  In  this  appendix,  the  initial 
subroutines  needed  to  define  Problem  (1)  with  Region  (1)  are  given 
for  the  Compressed  Jacobi  Conjugate  Gradient  Method. 


o n o o n on 


”T 


s 


s 


c »<•'<» 


c 

r 

C 


1H>9 


BEST  AVAIUBLE  COPY 
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ss  , 

,'A  M 

pppp  L 

EFFEI- 

PPRP 

s 

A 

A 

K*  J\  MM 

p p l_ 

F 

p 

A 

M M M 

P D 1 

r 

p 

SS'i 

A 

A 

^t  M M 

pppp  t 

FEF 

pppp 

A A A A u 

M M 

p 

f- 

p p 

A 

,Jt  M 

p 1. 

f. 

p p 

ss  , 

ft 

A 

Tv1  M 

p LLLLL 

FHEFi- 

p 

OR. 40. 

14  /4  J N|  77 

w 


p^Oi'Wi.M  r ir.GT'  r ( INiPijT  .OUTPUT  ) 


-’!-:  yl  GPIOX  (i'l  ) .GPIOY  (?1  ) .COPF  (^41  ,4)  .wO^'K 
-E-AL  RVALUM^) 

T iTP:GFP  GTYPF  (PI  .?!  ) .Wi^FO  (44  1 ) . IMV-rTA  (441  ) 

HA  T A TPT4/4L OUTPUT/. [POP/SL  TNPnT/ 

rORMON  W0P-<SP 

HEGIN:  COMMON  OECK  - ITPACK 

COMMON  / TTP'CK  / NC'’WPTG.NPOPTS.NPKPTS»NPOPP1  « 

A CMUF  t ';(.iUF  ,/PT  A ,rPSI . F ♦ OAmm  a . Pmq  ♦ S I OF  . 

R HAI.T.CA4FT  .CHANGE  » 

C PP.nELNNM.DELSNM.unNM.TFSTl  .Qa.(jT. 

0 IN.IS.ITmax, 

E I*-TR.TPieR, 

F OMEGA, fpfcTPA. PET ABAP.OmEGCHG 

HALT.CASFI .change .OmpgChg 
COMMON  OECK  - ITPACK 

COMMON  deck  - FLLPACK 


( 1 GP3)  .UNKNVN (44  1 ) 


LOGICAL 
END  : 

RFGIN: 

COMMON 

COMMON 

\ 

COMMON 

COMMON 

COMMON 

COMMON 

COMMON 


A 

p 

C 


I iTEGER 


logical 


A 

R 


REA( 


END 


/ PNOPM  / ipiecf.nhound.nrndpt 

/ CONSTS  / IPACM . IPACK?. TPACXH. INSIDE. HOP?. VERT. HOTh, 

corner. inter 

/ CONTHL  / DEoUG. LEVEL 
/ CPDF  / CUXX.CUXY.CUYY.CUX.CUY.CU 
/ F(iFORM  / NUM -iFQ.NUMCOE 
/ EONOEX  / NROW.NCOL 

/ PROP  / DIM?. PIM3. POISON. LAPLAC.CONSTC.SEL  FAD. CROSST, 
DIPICh.NFIIMAN.MIXEO.  AX.BX,  AY.RY.  AZ.RZ. 

NGR I DX , NGR 1 D Y . NGR I D Z . UN  I F R M . H X . HY , HZ . 
ELLP77.RFCT AN 
HORZ .VERT .POT H.CORNFR. 

P I E CE . HP  T YPE . RNP I GH , BGR I n 

niM?,OIM3.POISON.LAPLAC.CONSTC.SELFAD.CROSST.DIPICH. 
NFUMAN.MI XFn.UNIFRM,0F8UG.FLLP77.PECTAN, 

PTNAY 

AX.pX. AY. by. A/. RZ.HX.H Y. HZ. CUXX.CUXY.CUYY.CUX.CUY.CU. 

X ROUND. YBOUND. BP ARAM 
COMMON  DECK  - ELLPACK 


^EADdROP.TO)  ICASFS 
DO  GO  1 )KLM  = 1 . ICASFS 

INTEPFACE  i:  initial  situation 


n n o o n o o ('  o o o r>  r.  n o o o n n i*'-  o 
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s*-mplfp  - ? 


r 


[ 

[ 


l>EHU6  = .T»UF. 

LE'/FL  = ^ 

'JGSOXP  = p\ 

\G»OYD  = ?1 
-’XNCOE  = 4 
’MNHJ  = 44 1 

ITMAX  = son 

/f-:T4  = .000001 
= ,000001 
r^uE  = n.o 

CASE!  = .FAISP. 

K = .7>-. 

define  problem  on  OUTPUT  FIl.E 
WRITE ( IPTR.HO) 

WRITE ( IPTR»y0)  F.CMUE»2FTA»FPSI 

■IHHHHHHHHf-lHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHIHHHHHHHt 

interface  ?:  OO'-IATN  PROCESSING 


I 

t 

! 

I 


fr 


CALL  RFGIONIGTYPE .GRinx»NGROXO.GPlOY»NGROYD) 

I->)TEPFACE  3:  equation  GENERATION 

CALL  FIVEPT  (GTYPE»GRIDX,NGRDXD»GRinY,NGRDYD.COEF»*^XNCOE  » 
A f^XNEO) 

imterfacf  4:  equation  Indexing 

CALL  WROPD  {NDXEQ»MXNEQ«  INv/mDX) 


interface  s:  equation  -.olution 

CALL  INTUNK  (GT  Y^F  *GP  inx  ,NGRnXD»GRTOY»NGRDYD*MXNE{J» 

A HNKNWN) 

C ALL  C.tCGlGT  YPF  ♦GRinX,NGRnxn»ORlDY.NGPDYD«COEF»MXNCOE  »PXNF(J, 
A NOXFQ.i'NKNWN) 


C 

r 


P'TERFACE  b:  EiyTPUT 


NO  so  I I = I.NGRPTS 
IX  = MOD( I J-1 .NGWIUX ) ♦ 1 

JY  = ( IJ-I X) /NGRIOX  ♦ 1 
GO  TO  (10.20«30)  GTYPE(IX.JY) 
10  CONTINUE 


wQRKSPdJ)  = 
r,o  TO  40 


I 


TRUE (GRtOX ( IX ) .GRIOY ( JY)  ) 
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1 


BEST  AVAIUBIE'COPY 


?0  CONTINUE 

CALL  BCONO  ( inilMMY.GPinx  ( IX)  .GPIOY  ( JY)  .PVALUS) 

WOPKSP(TJ)  = RVALUS (4) /MYAl  US ( 1 ) 

GO  TO  40 

30  continue 

worksp(t  ))  = 0.0 

40  <~OMT  INUE 

WOPKSP ( T J*NCPPTS)  = WOPKSP(IJ)  - UNKNWN(IJ) 

SO  CONTINUF 

TPUNORm  = UTOV  (GTYPf  «NGRDXn.NGPOYntCOEF»MXNCOF  .^^XNEO  . NDXEO  » 

A WORKSP ( 1) »WORKSP ( 1 ) , 1 tNGRPTS) 

FRRNOR"  r UTOV  (GTYPE  *NGROXn,NGRDYO»COEF (MXNCOE  »mxnEO»NDXFQ« 

A WORKSP ( 1 ♦NGRRTS) .WORkSP ( 1 ♦NGRPTS) . 1 ♦NGR^TS) 

TRUNORM  = SORT (TRUNORM) 
i-RRNOR  -1  - SORT  (EPRNORN  ) 

WRITF  (IPTR^lOO)  TRUNORM. FRRNORM 

CALL  VOUT  (GRIOX»NGRDXO*GRir'Y.NGRnYD»WORK5P  ( 1 ) .NGRPTS) 

^0  CONTINUE 
70  format (IS) 

40  format  ( 1 HI  .////30X  ,*THF  method  hfimg  USED:  CJCG<>»/30X. 

A «the  ordering  heing  used:  REO-HLACKo. 

B /loX.^CASE  II  OF  THE  ADAPTIVE  PROCEDURE<> » //30X . 

C <>BOUNDARY  values  ARE  SET  TO:  0.0*. /iOX. 

0 » initial  SOLUTION  IS:  ().0<>./30X» 

E <^TEST  PRO^^l.EM  NO.  ?*) 

RO  format (///30X . oSTART ING  ITERATIVE  PARAMETERS  ARE:»»/3SX. 

A i>F  =».3X.F15.H»/3bX.<>CMUE 

B F1S.8./3SX»«-7ETA  =#  . E 1 S . B . / 3S  X . RS  I =*fF15.4) 

100  FORMAT ( IHI .//30X .«0  TO  1/2  NORM  OF  TRUE  SOlUTION  =* . E 1 S . H , / 30 X . 
A ^>0  TO  1/2  NORM  OF  TMF  ERROR  =<> . E 1 S . H . // ) 


FMO 


SUBROUTINF  PDF (X»Y»CValuS) 
real  CVALUS(7) 

DATA  PI/3. 141592653^8979/ 

C 

C TWO  OlMF^iSlONS 

C VALUES  OF  EOUATION  COFFFICSINTS  AT  (X.Y)  IN  hROER: 

C UXX.UXY.UYY.UX.UY.il. RIGHT  FIDE 

C 

FXY  = FXR(X<>Y) 

CVALUS  ( 1 ) = FXY 

CVALUSl?)  = 0.0 

CVALIJS(3)  = l.O/F.XY 

CVALUS(4)  = 0.0 

CVAI.US(5)  = 0.0 

CVALUS(6)  = -l./(l.  ♦ X ♦ Y) 

TSX  = STtPI^X) 

TSY  = SIm(PI<»Y) 

TCX  = COS(Plox) 

TCY  = C0S(PI»Y) 

E2XY  = EXV»FXY 


ODO  nooo  o ooooo 
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"1 


temp  s pi*  ( X*TSX*TCY  ♦ 3.*Y*Ei?XY*TCX*TSY) 

TEMPI  s TSX*TSY* ( (?.*Y*Y-PI*PI ) *E?XY  - PI*PI  -F X Y/ ( 1 , ♦X  ♦ Y ) ) 
CVAIUS(7)  = TFMP^TEMPl 

PPTURNJ 

FnO 


'SUBROUTINE  HCOND  ( I »X»Y»BVALUS> 
peal  BVALUS(c*) 

values  of  HOUNDAPY  conditions  COFFFICEINTS  AT  (X,Y) 
IN  the  OWDER: 

u«ux»uy*right  side 


BVALUS ( 1 ) 

= 1.0 

BVALUS (?) 

= 0.0 

BVALUS(3) 

= 0.0 

BVALUS(<») 

= TRUE.  (X*Y) 

RE  turn 

^ . r* 

FmH 


FUNCTION  ■^PXUNK(X»Y) 

initial  approximation  TO  UNKNOWN  VALUES 

APX'iNK  = 0.0 

Rt-  TURN 
ENO 


FU^!f'TION  TRUE(X«Y) 

data  pi /3,  U15'^^653S8979/ 

TRUE  SOI  UTION 

TRUE  = EXP{X*Y)*SIN(PI*X)*SIN(RI*Y) 

PkTURN 

EnO 


BEST  AVAIIABIE  COPY 


T 


^ 
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Appendix  3 

The  Subroutine  REGION 

REGION  Is  a subroutine  which  superimposes  a grid  of  size  h 
on  a region  defined  by  closed  contours.  This  routine  constructs 
a two-dimensional  Integer  array  over  the  smallest  rectangle  cir- 
cumscribing the  possibly  Irregular  region  and  denotes  each  grid 
point  with  Integer  values,  namely,  +1  for  Interior  points,  +2  for 
boundary  points,  +3  for  exterior  points.  To  utilize  REGION,  the 
vertices  defining  the  boundary  of  each  contour  In  the  particular 
region  are  specified  and  ordered  so  that  the  interior  of  the  region 
always  lies  on  the  left.  The  x and  y coordinates  of  the  end- 

i 

points  of  each  consecutive  line  segment  defining  a contour  are  given  j 

as  input  data.  The  permissible  line  segments  are  those  In  an  j 

arbitrarily  chosen  xy-plane  which  are  parallel  to  the  x-axls  or  the  | 

y-axls  or  which  form  a 45®  angle  with  an  axis  whose  endpoints  are  I 

grid  points  for  the  prescribed  h.  | 

While  REGION  was  originally  developed  several  years  ago,  it  | 

has  been  modified  and  Improved  recently.  This  recoding  has  removed  ^ 

restrictions  such  as  the  limits  on  the  number  of  allowable  vertices 
and  on  the  number  of  possible  contours.  REGION  is  now  coded  In 
standard  Fortran  with  an  improved  data  structure  and  with  optimized 
code  where  possible.  The  subroutine  REGION  is  now  compatible  with 
code  specifications  outlined  in  the  ELLPACK  Contributor's  Guide  [ 5] . 

Hence,  it  is  being  utilized  at  UT  Austin  as  an  ELLPACK  module  to 
perform  domain  processing.  While  REGION  is  somewhat  limited  with 
regard  to  the  types  of  domains  it  can  process,  it  does  work  success- 
fully on  very  compl icated regions  with  a number  of  "holes"  in  them-- 
all  defined  with  horizontal,  vertical,  and  45°  line  segments. 

As  an  illustrative  example  of  the  use  of  subroutine  REGION, 
consider  the  two-contour  region  (4). 
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(1.0) 


The  contours  are  defined  by  the  labeled  endpoints  of  each  line 
segment.  The  input  data  is  read  using  format  1615  and  consists  of 
the  number  of  contours,  the  number  of  vertices  for  a contour 
followedby  the  coordinates  of  the  vertices  from  their  rational  form, 
i.e.,  xl  x2  yl  y2  designate  vertex  (xl /x2 , y 1/y 2) . The  final  input 


data  is  the 

grid 

spacing  h 

in 

rational 

form. 

For 

example , 

if 

h = 1/20  is 

specified , 

then 

the 

input 

data 

would 

be 

as  follows 

• 

10  0 

10 

10 

10 

0 

10 

5 

10 

5 

10  5 

10 

10  10 

10 

10  2 

10 

2 

10 

4 

10 

4 

10 

4 

10  4 

10 

1] 

I' 


il 


The  printed  output  from  REGION  and  the  structure  of  the  Integer 
array  GTYPE  for  this  input  data  is  as  follows. 


BtSrAVMlABlt  CO?Y 
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i>  » O «•««««« ‘t> 

21 

^2^222^^???3333333333 

?0 

?0 

2111111111 P3333333333 

IW 

19 

211111111 1 ?33333i3333 

1^* 

1 H 

21  1 1 1 1 1 1 1 12333333333  i 

17 

1 7 

211111111 123333333333 

16 

211111111 123333333333 

IS 

15 

21111111 1123333333333 

14 

14 

21111111 1123333333333 

13 

13 

211111111 123333333333 

12 

12 

2111111111 23333333333 

1 1 

1 1 

21111111 1123333333333 

10 

10 

2111111111123333J3333 

* • • • *****  • • 

• * 

9 

211122222111233333333 

H 

• • * 

rt 

211123332111123333333 

7 

*^-«*  *•• 

• • • * 

7 

21 11233321 1 1 112333333 

*•••*  *•• 

. . . .« 

6 

211123332111111233333 

'■> 

* • • • *****  • • 

«• 

5 

211122222111111123333 

a 

t t t t * 

4 

21111111111111111233^ 

■ A ■ • t * 

3 

211111111111111111231 

^ 

2 

211111111111111111123 

1 

1 

222222222222222222222 

REGION  Printed 

Output 

Integer  Array  GTYPE 

The  printed  output  from  REGION  for  the  other  five  test  regions 
follows . 


BtSl  AVAllABlt  COP'^ 
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21 

«»««#«««««««»«««««'{»«« 

(2) 

21 

« o « o « 

20 

t 1 1 • T 1 1 - ^ 

20 

o 

^ * 

*...* 

1^^ 

14 

o 

^ » 

IH 

O - I - T “ - 

7 1 T - I - 1 

IH 

o ^ 

^ « 

^>  o 

• • • 

17 

. ^ j ^ 

17 

o 

O O’ 

• • • 

16 



- ^ . . . . ^ 

lb 

o ^ 

o o 

lb 

^ ^ ^ ^ 

lb 

o ^ 

• • • 

o o , ^ . 0- 

14 

^••TfttT-- 

. ^ ^ 

14 

o ^ 

• • • 

,o  ^••■•••^ 

13 

^41- 

13 

o ^ 

• • • 

ooooo  o 

12 

1 - 1 

12 

o 

• • • 

^ . . . . - . ^ o 

11 

- ...  O’ 

11 

o 

• • ♦ 

. o 

10 



1 ‘11*^^ 

10 

0 ^ 

• • • 

o 

9 

^1  t 1 1 1 f “ * 

O’ 

9 

• • • 

o 

U 

•0 

O’ 

4 

* • 

• • • 

o 

7 

^ . , - , - , o 

7 

• • • 

b 

...  o 

b 

o ^ 

• • • 

S 

^ ^ ^ ^ o 

b 

o ^ 

• • • 

4 

*liTttt--“ 

- 1 T 1 ^ 

4 

o ^ 

• • • 

ff“r--r  i*ri 

J 

o 

3 

-ci-  ^ 

• • • 

i--i  * 

2 

2 

<}•  ^ 

• • • 

« ’tnnf  iHi- ii' in^ 


2\  * 
20  » 
1^  ^ 
IH  o 
17  ^ 
16 

lb  <> 


«««««»««» 


••••••••^ 

••••••••^ 

••••••••^ 

• •••••••'^ 

••••••••^ 


»«««»««««««««««»«•«»« 
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(6) 


H = 1/40 


37 

36 

35 

34 

33 

3^ 

31 

30 

26 

27 

26 

25 

24 

23 

2/ 

21 

20 

IH 

17 

16 

15 

14 

I j 
12 

I I 
10 

V 

« 

7 

6 

5 

4 

3 

2 

I 


* •••»••••••********'""“*** 




•innnnnnnnnnnn>ejnn>,,,,,,,,,, 



* 

.-.O 


iKHi  • • • • 

. . .« 

• * 

* • 

• 

. . .* 

■ft  ^ 

. . . 

. . .* 

^ • • • • • • 

* • 

^ ^ . 

...* 

• • • . . 

^ t T T * * - 

. . .* 

« 

• * 

* • 

^ ^ - 

• • * 

^ ^ 

• * 

««« 

O - - ^ - 

• 

* • • 

^ «■ 

t t t t « * t ^ 

• • 

• • ♦ • 

• • 

• • • 

^ ^ 

<> 

• • 

• • • • 

• • 

• • • 

1 t t t t ^ 

• • • • 

• • 

• • • 

• • • • 

• • 

• • • 

. . .* 

^ • • • 

• • 

• • • 

• • 

• • • 

• ^ 

«•  ^ 

• • 

• • • 

• • 

O 

For  completeness,  we  now  give  the  listing  of  subroutine  REGION 
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RFGION 


1 


RRRR 

EEEEE 

GGGG 

III 

000 

N 

N 

R R 

E 

G 

I 

0 0 

NN 

N 

R R 

E 

G 

I 

0 00 

N 

N N 

RRRR 

EEE 

G GG 

I 

0 0 0 

N 

NN 

R R 

E 

G G 

I 

00  0 

N 

N 

R R 

F 

G G 

I 

0 0 

N 

N 

R R 

EEEEE 

GGGG 

III 

000 

N 

N 

11 

.42.52 

07  JUN 

77 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


«^UHR0UTTNE  RFGION  ( GT  YPE  tGR  lOX  .NGROXD  * GR  lOV  » NGROYO ) 


FUNCTION  : SUPERIMPOSES  A MESH  OF  SIZE  HX=HY=H1/H2  ON  A REGION 
DEFINED  RY  CLOSED  CONTOURS  AND  CONSTRUCTS  AN  INTEGER 
ARRAY  WHICH  DESCRIBES  EACH  MESH  POINT  ON  THE  SMALLEST 
rectangle  circumscribing  the  possibly  irregular  REGION 
AS  AN  INTERIOR  POINT  (♦!).  AN  BOliNOARY  POINT  U2)»  OR 
AN  EXTERIOR  POINT 


USAGE  : CALL  region  (6TYPE*GRIDX»NGR0XD»GRI0Y»NGR0YD) 

PARAMETERS  : 

GTYPE  - GTYPE  IS  AN  NGROXD  BY  NGROYD  INTEGER  ARRAY  USED  TO 

INDICATE  THE  Type  OF  POINT  ON  THE  GRID.  THE  NUMBERS 
!♦  2«  OR  3 INDICATE  RESPECTIVELY  INTERIOR*  BOUNDARY* 
OR  EXTERIOR  POINTS  OF  The  GRID. 

NGRDXD  - N6RDXD  IS  THE  ROW  DIMENSION  OF  THE  ARRAY  GTYPE  AS 
SPECIFIED  IN  THE  CALLING  PROGRAM. 

NGRDYD  - NGROYD  IS  THE  COLUMN  DIMENSION  OF  THE  ARRAY  GTYPE 
SPECIFIED  IN  THE  CALLING  PROGRAM. 

GRIDX  - GRIDX  IS  AN  ARRAY  OF  LENGTH  NGRDXD  DIMENSIONED  IN  THE 
CALLING  program.  UPON  LEAVING  REGION  IT  CONTAINS  THE 
X COORDINATES  OF  THE  MESH  LINES  STARTING  IN  THE 
LOWER  LEFT  HAND  CORNER. 

GRIDY  - GRIOY  is  AN  ARRAY  OF  LENGTH  NGROYD  DIMENSIONED  IN  THE 
calling  program.  upon  leaving  region  it  CONTAINS  THE 
Y COORDINATES  OF  THE  MESH  LINES  STARTING  IN  THE 
LOWER  LEFT  HAND  CORNER. 


OTHER  PARAMETERS  PASSED  IN  LABELED  COMMON  ARE*. 


level 


0 NO  PRINTING  FROM  REGION 

1 THE  INPUT  DATA  ONLY  IS  PRINTED. 

? THE  GRAPH  OF  THE  REGION  ONLY  IS  PRINTED. 

3 PRINT  BOTH  input  DATA  AND  GRAPH  OF  REGION 


DEBUG  IS  A LOGICAL  DEBUGGING  PARAMETER.  IF  TRUE  THEN  LEVEL 
IS  RESET  TO  3.  IF  FALSE  NO  ACTION  IS  TAKEN. 


BEST  AVAIUBLE  COPr 
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PFGION 


p 


C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


BEST  AVAIWBlE'COPy 


NGPIOX  IS  THE  NUMHER  OF  MESH  POINTS  IN  THE  X-OIRECTION  OF  THE 
CIRCUMSCRIBED  RECTANGLE.  THIS  IS  COMPUTED  IN  REGION. 


NGRIOY  IS  The  number  of  mesh  points  in  the  Y-DIRECTION  OF  THE 
CIRCUMSCRIBED  RECTANGLE.  THIS  IS  COMPUTED  IN  REGION. 


NGRPTS  is  the  number  of  TOTAL  MESH  POINTS  OF  THE  CIRCUMSCRIBED 
rectangle.  THIS  IS  COMPUTED  IN  REGION. 

HX.HY  ARE  THE  MESH  SIZE  FOP  THE  GRID.  REGION  READ*:  h1  AND 
H2  FROM  DATA  AND  COMPUTES  HX=HY=H1/H2.  HX  AND  HY  ARE 
THEN  RETURNED  TO  THE  CALLING  PROGRAM, 

AX»BX  ARE  the  minimum  AND  MAXIMUM  VALUES  OF  THE  X COORDINATE. 

REGION  COMPUTES  THESE  AND  RETURNS  THEM  IN  LABELED  COMMON. 

AY. BY  ARE  THE  MINIMUM  AND  MAXIMUM  VALUES  OF  THE  Y COORDINATE. 

REGION  COMPUTES  THESE  AND  RETURNS  THEM  IN  LABELED  COMMON. 


PARAMETER‘S  USING  BLANK  COMMON  ARE: 
L IS  DESCRIBED  BELOW 


other  parameters  : 

H1.H2  ARE  parameters  INDICATING  THE  UNIFORM  MESH  SIZE 

IN  rational  FORM.  THESE  AS  WELL  AS  L ARE  READ  IN  AS 
DATA  AND  ARE  DESCRIBED  FURTHER  BELOW. 


THE  DEFINITION  OF  THE  REGION  IS  AS  FOLLOWS: 

KN  = number  OF  CONTOURS  IN  THE  REGION. 

L(5.K)  = NUMBER  OF  VERTICES  ON  THE  K — Th  CONTOUR  --  1 OR  MORE 
L(1.I)/L(2.I)  = THE  X COORDINATE  OF  THE  I-TH  VERTEX. 

L (3. 1 ) /L  (*..  I ) = THE  Y COORDINATE  OF  THE  I-TH  VERTEX. 

H1/H2  = THE  mesh  SIZE  TO  BE  CONSIDERED. 

notes:  (1)  THE  ARRAY  L USES  BLANK  COMMON  FURNISHED  IN  THE 
CALLING  PROGRAM 

(2>  THE  INPUT  DATA  FOR  THE  VERTICES  MUST  BE  ORDERED  SO 
THE  INTERIOR  OF  THE  REGION  ALWAYS  LIES  TO  THE  LEFT 
(3)  REQUIRED  SUBROUTINES  --  RTREG. I ABS.MOOS 


WRITTEN  OR  MODIFIED  BY 


DATE 


ROGER  G. 
ROGER  G. 
JAMES  0. 
ALKIS  J. 


grimes  and  DAVID  R.  KINCAID 
grimes  AND  DAVID  R.  KINCAID 
SULLIVAN  AND  DAVID  R.  KINCAID 
MOURADOGLOU  and  JOHN  H.  DAUWALDER 


JUNE  I*)?? 
FEBRUARY  1977 
FEBRUARY  197ft 
APRIL  1967 


CENTER  FOR  NUMERICAL  ANALYSIS/COMPUTaTION  CENTER 
university  op  TEXAS  AT  AUSTIN 


•>on  ooooooo  ooooonnooooooo 


references: 


BESl  AVAILABLE  COPY 


81 

REGION  - 3 


(1)  KINCAIDtDAVID  R.  AND  ROGER  G.  GRIMES,  *****  CNA  REPORT  *** 


(2)  ALKIS  J.  MOURADOGLOU  AND  JOHN  H.  OAUWALOER*  rREGION 

SPECIFICATION  ROUTINEE,  V2  UTEX  REGION,  UT V2-0 1 -CC026 , 
COMPUTATION  CENTER,  UT-AU'=TIN,  APRIL  19G7, 


(3)  A.  J,  MOURADOGLOU,  ENUMEPICAL  STUDIES  ON  THE  CONVERGENCE 

OF  THE  PEACEMAN-RACHFORD  ALTERNATING  DIRECTION  IMPLICIT 
METHOOE,  TNN-67,  COMPUTATION  CENTER*  UT-AUSTIN,  JUNE  1967. 


real  GRIDX (NGROXD) *GRIDY (NGRDYD) 

INTEGER  GTYPE (NGRDX0,NGRDY0» 

INTEGER  R,S,RS1 ,SSAVE,SET,P,Q,FlAG,H1 ,H2 
COMMON  letter (3) ,NUMBER(3) ,L(5,1) 

***  begin:  common  deck  - ITPACK 
**•  END  : COMMON  DECK  - ITPACK 


***  BEGIN:  COMMON  DECK  - ELLPACK 
***  END  : COMMON  DECK  - ELLPACK 

number ( 1 ) =2 

NUMBER(2) =1 
NUMBER (3) =3 
LETTER ( 1 ) =1H* 

LETTER(2)=1h. 

LETTER(3)=lH 

MPl=NGROXn-l 

MQ1=N6RDY0-1 

READ  (IRDR»610)  KN 

READ  (IRDR*610>  KCOl 

L (5, 1 ) =KCOl 

READ  (IRDR*610)  <L ( 1 , J) ,L ( 2, J ) tL ( 3, J) , L ( ^* J ) , J= 1 ♦ KCOl ) 
IF  (KN.LE.l)  GO  TO  20 
00  10  K=2,KN 

READ  (IROR,610)  L(5,K) 

KC02=KC01*1 
KC01=KC01*L (5,K) 


***  IT  IS  ASSUMED  THAT  L(5,K)  > 0 FOR  K=1,2*...,KN  AND  THAT  KN  > 0. 

10  READ  (IRDR»610)  (L ( 1 » J) »L (2, J) »L (3* J) tL (4, J) , J=KC02,KC01 ) 

r 

C***  THE  OEAOING  OF  THE  COORDINATES  OF  THE  VERTICES  IS  WITH  1615  FORMAT 
C***  L«1*I>*L(?*I)*L(3,I),L(4,I),L(1»I*1)*L(2,I*1),L(3,I*1),... 

20  READ  (IROR*610)  H1,h2 
K J2a0 

IF  (DEBUG)  LEVEL=3 
IF  (M00(LEVEL,2) .EQ.O)  GO  TO  40 


• r*  T iNPtjT  DATA 


BEST  AVAILASIE  COPY  ■ “ 

WRITF  (IPTR.620) 

KJ2  = 0 

no  30  K=l,KN 
KJ1=KJ2*1 
KJ2=KJ2*L (5»K) 

WRITF  (IPTR»630)  K 
WRITF  (IPTR»640) 

WRITF  (IPTR»660)  (L ( 1 » I ) ♦L (2» I ) » I=KJ1 ♦KJ2) 

WRITF  (IPTR,650) 

30  WRITF  (IPTR»660)  (L (3» I ) *L (4* I ) ♦ I=KJ1 »KJ2) 

K J2=0 

40  00  60  1=1, KN 
KJ1=KJ2^2 
KJ2=KJ2*L (5»  I ) 

IF  (KJ1,GT.KJ2)  STOP  1 
DO  50  KS=KJ1,KJ2 
KS1=KS-1 

IF  (L(1,KS)*L(2,KS1) ,FQ.L(2»KS)*L(1»KS1) .0R.L(3,KS)*L(4,KS1) 

1 .FQ.L (4,KS) *L (3,KS1 ) .OR. I ARS ( (L ( 1 *KS ) «L (2,KS 1 ) -L ( I tKS 1 ) *L  <2, 

2 KS))*L(4,KS)*L(4tKSn)  .EQ.  IABS(L  (2»KS)«L(2,KSl)*(L(3,KS)*L(4 

3 ,KS1 )-L(3»KSl)«L(4,KS) ) ) ) ^0  TO  50 
WRITE  (IPTR,670)  KS,I 


STOP  2 

50  CONTINUE 
KJ1=KJ1-1 
DO  60  J=1.KN 

IF  (L(l»KJl)«L(2,KJ2) .FO.L ( 1 »KJ2) *L (2,KJ1 ) .OR.L (3»KJl ) *L (4,KJ2) 

1 .EQ.L (3,KJ2)*L {4,KJ1 ) .OR. IASS ( (L ( 1 »KJ1 )«L (2*KJ2)-L ( 1 ,K J2) *L (2,K 

2 J1  n«L  (4,KJ1  I4,kJ2)  ) .EQ.  I ABS  (L  (2, KJl ) *L  C2»KJ2)*  JL  (3»KJ1  ) *L  (4 

3 ,KJ2)-L(3»KJ2)*L(4,KJ1) ) ) ) GO  TO  60 
WRITF  (IPTR,670)  KJl.J 

<;T0P  3 
60  CONTINUE 


C«** 

C*** 

c**» 

c*** 

c*** 

c*** 

c 


PICK  THE  MAX  AND  MIN  OF 
ALL  CONTOURS. 

MINX1/MINX2  = MIN  OF 
MINY1/MINY2  = MlN  OF 
MAXX1/MAXX2  = MAX  OF 
MAXY1/MAXY2  = MAX  OF 


the  coordinates  of  THE  VERTICES  OF 

X COORDINATES 

Y COORDINATES 
X COORDINATES 

Y COORDINATES 


MAXX 1=L  (1,1) 
MINX 1=MAXX1 
MAXX2=L (2,1) 
MINX2=MAXX2 
MAXY1=L (3, 1 ) 


MINY1=MAxy1 
MAXY2=L (4,1 ) 
MINY2=MAXY2 
C 


c*** 

MINXl 

= MAXXl 

= X 

( 

1 , 

1 

) 

c*** 

MINX2 

= MAXX2 

= X 

( 

2 , 

1 

) 

c*** 

MINYl 

= MAXYl 

= X 

( 

3 , 

1 

) 

c*** 

MINY2 

= MAXY2 

= X 

( 

4 , 

1 

) 

C 


IF  (KCOl.LE.l)  STOP  4 
DO  100  I=2,KC01 
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70 


80 


90 


100 


110 


C 

c**» 

c 

120 

C 

c*** 

c 


130 


140 


150 

C 

c 


160 

170 


IF  (M4XX1*L(2.I) .GF.MAXX?*L(1»I) ) GO  TO  70 
MAXX1=L ( 1 ♦ I ) 

MAXX2=L (2. I ) 
r,0  TO  80 

IF  (MINX1*L  (2*  I ) ,LF.MINX2<»L  ( 1 » I ) ) GO  TO  80 
MINX1=l (1  .1  ) 

MTNX2=L (2.1 ) 

IF  (MAXY1*L (4. I ) .GE.MAXY2*L (3. I ) ) GO  TO  90 
MAXY1=L(3.I) 

MAXY2=L (4. I ) 

GO  TO  100 

IF  (MINY1*L (4. I ) .LF.MINY2*L (3. I ) ) GO  TO  100 
MINY1=l (3.1) 

MINY2=L  (4.1) 

CONT INUE 

IF  (H2*(MAXX1*MINX2-MAXX2*MINX1) ,LE.MP1*H1*MINX2*MAXX2)  GO  TO  110 
WRITt^  (IPTR.680) 

STOP  5 

IF  (H2*(MAXY1*MINY2-MAXY2*MINY1) .LE.MQ1»H1*MINY2*MAXY2)  GO  TO  120 
WRITE  (IPTR.690) 

STOP  b 

CHECK  That  all  boundary  points  are  integral  multiples  of  H 
DO  150  I=1.KC01 

AT  THIS  POINT  KCOl  = THE  TOTAL  NUMBER  OF  VERTICES. 

MM=H2* (L ( 1 . I ) *MINX2-L (2.1 )*MINX1 ) 

NN=h1«MINX2«L  (2.1) 
kk=mm/nn 

IF  (KK*NN.FQ.MM)  GO  TO  140 
WRITF  (IPTR.700) 

STOP  7 
L ( 1 . 7 ) =KK 

MM  = H2*  (L  (3.  I )*MINY2-L  ( 4 . 1 ) <»M  I NY  1 ) 

NN=H1*MINY2*L (4. I ) 

KK=MM/NN 

IF  (KK»NN.NE.MM)  go  to  130 
L (2. I ) =KK 
CONTINUE 


determine  ngriox  and  NGPIDY, 


njGRXMUl  (1.1) 

NGRYM1=l (2.1 ) 

IF  (KCOI.LT.2)  STOP  10 
DO  170  I=2.KC01 

IF  (NGRXMI .GT.L ( 1 . I ) ) GO 
NGRXM1=L(1.I) 

IF  (NGRYMl .GT.L (2. I > ) GO 
NGRYM1=L (2. 1 ) 

CONT INUE 
NGRI0X=NGRXM1 .1 
NGPXP1=NGRI0X.1 
NGRI0Y=NGPYM1 .1 
NGRYPUNGPlOY.l 


TO  160 
TO  170 

Btsi  AVWlABlt  C0P'( 
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BEST  AVAIUBIE'COPY 

c 

c***  SET  THE  ARRAY  iGTYPEE  Tn  ZERO. 

C 

00  180  J=1»NGRI0Y 
DO  IHO  I=1»NGRIDX 

180  GTYPE(I,J)=0 

C 

c***  define  THE  BOUNDARY  POINTS  WHICH  ARE  NOT  VERTICES 

C***  IN  THE  ARRAY  GTYPE 

C 

KJ2  = 0 

00  480  K=1»KN 
KJ1=KJ2^1 
KJ2=KJ2^L (5»K) 
no  480  J=KJ1»KJ2 

IF  (J.NE.KJl)  GO  TO  190 
IPP=1 ♦L ( 1 .KJ2) 

TQP=1*L(2«KJ2) 

GO  TO  200 

IQO  IPP=1  1 • J-1 ) 

IQP=1»L(2»J-1) 

200  IP=1*L(1»J) 

IQ=UL(a*J) 

IF  (J.NE.KJ2)  GO  TO  210 
TPN=1*L(1»KJ1) 

IQN=l4.L(2*KJl) 

GO  TO  220 

210  !PN=l4L(l.J*l) 

I0N=1*L(2»J*1) 

220  CALL  RTREG  ( IPP » I OP » IP » 10 »R ) 

CALL  RTREG  ( IPN» IQN» IP, IO»S) 

IF  (GTYPE(IP,I0)  .EQ.O)  GO  TO  230 
WRITf  (IPTR,710)  J,K 
STOP  11 

230  M0DR=M0D(R*4,8) 

M00S=M0D(5,4,8) 

IF  (R.LE.S)  GO  TO  2b0 
IF  (MODR.LF.MODS)  go  to  240 
GTYPEdP,  IQ)=,1 
GO  TO  280 

240  GTYPE(IP,IO)=,10 

SO  TO  280 

250  IF  (R.NE.S)  GO  TO  260 

wRITr  {IPTR,720)  J,K 
STOP  12 

260  IF  (MOOR. LT. MODS)  GO  TO  270 

GTYPf{IP,IQ)=-1 
GO  TO  280 

270  GTYPf ( IP, 10) =-10 

280  IF  (J.NE.KJl)  GO  TO  290 

RSI=R 
SSAVf=S 
GO  TO  480 

290  IF  (lARS(IPP-IP) .LE.I.ANO.IABS(IQP-IO) .LE.l)  GO  TO  470 

MOOR=MOD (R*4,0) 

M0DS=M0D(SSAVE*4,8) 

IF  (R.LE.SSAVE)  go  to  300 


o o n o o o 


tBI  AVAIUBIE  COPY 

SET=^10 

IF  (MOnR.GT.MODS)  SET=*1 
GO  TO  310 
300  SET=-10 

IF  (MOOR. GT. MOOS)  SET=-1 
310  IPM1=IaBS(IPP-IP)-1 

IQMUIaBS  ( IQP-IQ)-1 
IF  (IPP.GE.IP)  GO  TO  370 
IF  (lOP.GF.IQ)  GO  TO  330 
00  320  N=l,IPMl 
320  GTYPF ( IPP^N. IQP>N) =SET 

GO  TO  470 

330  IF  (lOP.NE.IQ)  GO  TO  350 

00  340  N=l,IPMl 
340  GTYPE(IPP'»N.IQ)=SET 

GO  TO  470 

350  00  360  N=1.IPM1 

360  GTYPF ( IPP.N. IQP-N) =SET 

GO  TO  470 

370  IF  (IPP.NF.IP)  GO  TO  410 

IF  (lOP.GF.IQ)  GO  TO  390 
no  340  N=l«IOMl 
340  GTYPE ( IP* IOP*N) =SET 

GO  TO  470 

390  00  400  N=l.IQMl 

400  GTYPF ( IP* IQP-N) =SET 

GO  TO  470 

410  IF  (lOP.GF.IQ)  GO  TO  430 

00  420  N=l,TPMl 
420  GTYPE(IPP-N»I0P*N)=SET 

GO  TO  470 

430  IF  (lOP.NF.IQ)  GO  TO  450 

00  440  N=l . IPMl 
440  GTYPF ( IPP-N* 10) =SET 

GO  TO  470 

450  00  460  N=l,IPMl 

460  GTYPE(IPP-N*IQP-N)=SET 

470  SSAVE=S 

IF  (J.NE.KJ2)  GO  TO  480 

IPP=IP 

IQP=IO 

IP=IPM 

IQ=IQM 

P=RS1 

J = J*1 


***  THE  DO  LOOP  INDEX  J HAS  JUST  BEEN  REDEFINED  INSIDE  THE  LOOP. 


ooo  ooo  ooo 
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BEST  AVAILABLE  COPY 

IF  (GTYPE ( I ♦ J) .EQ.O)  GO  TO  490 
F|.AG=1 
GO  TO  510 

490  IF  (Flag. EQ.O)  go  to  500 

Fl.AG  = 0 
SFT  = 2 

IF  (GTYPE (I-l .J) .GE.O)  5ET=3 
500  GTYPE(I»J)=SET 

510  continue 

**  CHECK  CONSISTENCY  OF  CONTOUR  ORIENTATIONS  SCANNING  RIGHT  TO  LEFT. 

FLAGrO 
SFT  = 3 

DO  540  K=1*NGRIDX 

i=ngrxpi-k 

IF  (GTYPEdtJ)  .EQ.2)  GO  TO  520 
IF  (GTYPE ( I tJ) .EQ.3)  GO  TO  520 
FLAG=1 
GO  TO  5A0 

520  IF  (FLAG. EQ.O)  GO  TO  530 

FLAG=0 
SET=? 

IF  (IA8S(GTYPE(I*1»J) ) .EO.l)  SET=3 

530  IF  (GTYPE( I« J) .EQ.SET)  GO  TO  540 

WRITF-  (IPTR»730)  I.J 
STOP  13 
5A0  CONTINUF 


**  SET  THE  VALUES  OF  THE  ARRAY  EGTYPEr. 


DO  550  J=1*NGRIDY 
DO  550  I=ltNGRIOX 
K=GTYPF ( I ♦ J) 

IF  ( (K.NE.2) .AND. (K.NE.3) ) K=1 
550  GTYPF(I,J)=LETTER(K) 

IF  (LEVEL. LT. 2)  GO  TO  570 

PRINT  THE  REGION  IDENTIFICATION  GRID. 

WRITE  (IOTR»620) 

WRITE  (IPTR»740)  H1.H2 
ME  tA0=NGRlDX 
DO  560  K=1»NGRI0Y 
j=NGRYP1-K 

560  WRITE  (IPTR»750)  J« (GTYPE ( I ♦ J) ♦ 1=1 *MAXAB) 

570  CONTINUE 

DO  580  J=l»NGRIOY 
DO  580  I=ltNGRIOX 

IF(GTYPE(I»J).EQ.LETTER(1))  GTYPE(I»J)  =NUMBER(1) 
IF(GTYPE(I»J) .EQ.LETTER(2) ) GTYPE(I»J)  = NUMBER(2) 
IF(GTYPE(ItJ)  .EO.LETTER(3)  ) GTYPEdtJ)  =NUMBER(3) 
580  CONTINUE 

ngrpts=ngriox*ngrioy 
HX=FLOAT (Hi ) /FLOAT (H2) 
hy=hx 

AX=FLOAT (MINXI ) /FLOAT (MINX2) 


o n o n o o 


BEST  AVAIUBLE  COPY 


87 

PFGION  - 9 


BXsFLOAT (MAXXl ) /FLOAT (MAXX2) 

AY=FLOAT (MINT  1 ) /FLOAT (MINY2) 

BY=FLOAT (MAXYl )/FLOAT (MAXY2) 

DO  590  ISET=1 .NGRIDX 
) GRIOX ( I SET) = AX ♦FLOAT ( lSET-1) *HX 
DO  600  ISET=1 ,NGRIDY 
) GRIDY(ISET)=AY^FL0AT(ISFT-1)*HY 
RETURN 

) format  (1615) 

) format  (1h1«////) 

) FORMAT  (///25Xtl3HCONTOUR  NO.  tlSt/) 

I FORMAT  (/»1H0»AX.25HX  COORDINATES  OF  VERTICES) 

) FORMAT  (/f 1H0*4X*25HY  COORDINATES  OF  VERTICES) 

) FORMAT  (lOXt8(lX»I4»lH/»I4»2X) ) 

) format  (1H2»///»10X»23HTHE  COORDINATES  OF  THE  tlS.lBH  TH  VERTEX  ON 

1 the  *I5tl2H  TH  CONTOUR  ♦/♦ I 0X»55HDEFINES  WITH  THE  PREVIOUS  VERTEX 

2 A segment  WHICH  MAKES  » / ♦ I OX ♦ 68HW ITH  THE  X-AXIS  AN  ANGLE  OTHER  TH 
3AN  N*(Pl/4)  WHERE  N IS  AN  INTEGER  ♦//) 

) FORMAT  ( lH2t///»10X»36HTOO  MANY  MESH  POINTS  IN  X DIRECTION  ♦//) 

) format  ( 1h2»///» IOX»36HTOO  MANY  MESH  POINTS  IN  Y DIRECTION  ♦//) 

) format  ( 1h2»///« 10X»20HH  IS  NOT  ACCEPTABLE  .///) 

) format  (1H2*10X»4HTHE  flStlBH  TH  VERTEX  OF  THE  tI5»24H  TH  CONTOUR 

ialready  set  »///) 

) format  (lH2tlOX*4HTHE  tlStlBH  TH  VERTEX  OF  THE  »I5*24H  TH  CONTOUR 
IGIVES  ^ = S ♦///) 

) format  (lH2«///»10Xt56HINCONSISTENT  CONTOUR  ORIENTATION  FOR  MESH  P 
lOINT  WITH  P =I3.2X*3HQ  =13) 

) format  (25X»31HTHE  REGION  IDENTIFICATION  GRID  ♦//♦30X*22H.  ARE  INT 
IERIOR  POINTS  ♦/♦30X»22H«  ARE  BOUNDARY  POINTS  ♦/* 30X , 28H ( BLANK ) ARE 
2 EXTERIOR  POINTS  ♦//♦35X»4HH  = . I4» lH/» I4»//) 

) format  (3X»I5»2X»101A1) 


subroutine  RTREG  (PI tOl tP2»Q2»RS) 
integer  pi  ♦P2.QI »Q2»RS 

IF  (P1-P2)  10»50«90 

»•»  PI  < P2, 

10  IP  (01-02)  20»30*40 
20  RS=5 
RETURN 
30  RS=4 
return 
40  oS=3 

return 


• PI  = P2, 


50  IF  (01-02)  60.70«80 


60  »S=6 

return 

70  WRITE  (iyTR«130) 

«;top  u 

PO  RS=2 
RETURN 
C 

C*»*  PI  > P2. 

C 


90  IF  (01-Q2)  100*110*120 
100  RS=7 
RETURN 

no  t>S=0 

RETURN 
120  RS=1 

return 

130  FORNAT  (1h2*10X*  42H(RTREG)  TWO  VERTICES  HAVE  THE  SAME  P AND  Q*//) 


C 


r 


A new  version  of  ITPACK  (August  1977)  has  added  capabili- 
ties to  the  version  covered  in  this  report.  These  capabilities 
are  a new  solution  method.  Symmetric  SOR  Partially  Adaptive 
(hereafter  referred  to  as  SSOR-PA) , a constant  coefficient  switch, 
and  an  adap t ive/ nonadap t ive  switch. 

SSOR-PA  is  similar  to  symmetric  SOR  Semi-iterative  (SSOR-SI) 
except  it  applies  the  adaptive  process  only  to  the  spectral 
radius,  SPECR.  To  use  SSOR-PA  a good  choice  of  CME  and  OMEGA  must 
be  known,  a priori. 

The  constant  coefficient  switch  is  a logical,  variable 
CONSTC.  If  the  user  wants  to  solve  an  Elliptic  Partial  Differen- 
tial Equation  with  constant  coefficients,  then  CONSTC  should  be 
set  to  .TRUE,  in  the  main  program.  This  will  allow  the  user  to 
set  MXNCOE*!  and  reduce  the  storage  allocation  by  3 full  vectors. 
The  array  COEF  would  then  contain  only  the  right-hand  side  of  the 
equation  and  the  constant  coefficients  would  be  saved  in  a labeled 
common  block. 

The  adaptive/non-adaptlve  switch  is  the  logical  variable, 
ADAPT.  If  the  user  already  has  a good  parameter  selection  and  does 
not  wish  to  change  parameters  adaptively,  ADAPT  should  be  set  to 
.FALSE.  This  switch  is  available  in  all  solution  methods  except 
SSOR-PA.  As  SSOR-SI  and  SSOR-PA  are  identical  in  the  non-adaptlve 
case,  this  option  was  added  only  to  SSOR-SI  as  it  required  less 
storage  allocation  than  SSOR-PA. 

Another  change  in  the  new  version  of  ITPACK  is  the  initial- 
ization of  scalars,  parameters,  and  switches  that  are  needed  in 
ITPACK.  To  Interface  with  the  ELLPACK  Control  Program  these 
variables  could  not  be  initialized  in  the  main  program  and  must 
Instead  be  initialized  from  an  input  file.  The  following  variables 
are  still  initialized  in  the  main  program: 


CONSTC 

DEBUG 

LEVEL 


90 


NGRDXD  CONSTC 

NGRDYD  DEBUG 

MXNCOE  LEVEL 

MXNEQ 

The  variables  which  are  initialized  from  the  input  file  are; 


ITMAX 

ZETA 

EPSI 

ADAPT 

CASEI 


F 

CME 

SHE 

OMEGA 

SPECR 


The  input  file  for  these  variables  consists  of  two  lines 
and  must  be  added  to  the  end  of  the  input  file  for  the  REGION 
subprogram  (see  Appendix  3).  The  first  line  is  the  same  for  all 
solution  methods.  ITMAX,  ZETA,  EPSI,  ADAPT  and  CASEI  are  read  in 
a 110,  2F10.2,  2L10  format. 

The  variables  on  the  second  line  of  the  input  data  line 
depend  on  the  solution  method  but  all  are  read  in  4F10.2  format. 
For  J-SI  or  RS-SI  the  variables,  in  the  order  they  appear,  are  F, 
CME,  and  SME.  For  SSOR-CG,  SSOR-PA,  and  SSOR-SI  the  variables  are 


F,  CME,  OMEGA,  and  SPECR.  In 

the  adaptive 

case 

(ADAPT-. TRUE.  ) 

RS-CG  and  CJ-CG  need 

no  second 

line,  while 

in  the  non-adaptlve 

case  (ADAPT-. FALSE.  ) 

CME  is  read  from  the 

second 

data  line. 

A sample  of  the 

above  for 

the  SSOR-PA 

solution  method  follows. 

Columns  - 10 

20 

30 

40 

50 

100 

. 000001 

.000001 

TRUE 

FALSE 

.75 

.95 

1.64 

o 

o 

ITPACK  would  then  set 

: 

ITMAX- 

100 

II 

i-n 

ZETA-. 

000001 

CME-. 95 

EPSI-. 

000001 

OMEGA-1 . 64 

ADAPT- 

.TRUE. 

SPECR-0. 0 

CASEI” . FALSE. 
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